format: text

text:
# 1 Introduction

## 1 Introduction

Financial markets are the object of considerable interest to an increasing number of physicists, because of their complex nature and the applicability of stochastic techniques to such dynamic systems. This is a relatively new area of interdisciplinary physics research, and is catching the attention of quite a few theoretical and condensed matter physicists, especially those working in the fields of complexity and chaos. The relationship between physics and finance is not something that started in recent years. In fact, as far back as 1900, Louis Bachelier proposed the random walk model of the stock market. Historically, much of the initial work on Brownian motion and fractals originated from studies of stock market behaviour, before being absorbed by mainstream physics. In the 1960s and 1970s these ideas became very popular and eventually led to the famous Black-Scholes option pricing formula.

Recently, many papers have appeared in which markets are treated as far-from-equilibrium dynamical systems. Econophysics is a quickly growing field of research, and work is being done on the scaling behaviour of exchange rates, “log-periodic” oscillations as crash precursors, dynamics of the interest rate curve, and market fluctuations, to name but a few examples. Some list of references (very incomplete) can be found in . For more detailed discussion on physical approaches to financial markets and economy we refer the reader to recent publications , ,  which provide a lot of references.

In this paper we present results of numerical simulations of two models - Minimal and Two-liquid, which had been formulated by two of us (D.E. and I.I.K.) in the paper on market microstructure . In our approach the market microstructure is described by a diffusion-annihilation model. Different states of market are described by different possible state of a system in which diffusion controlled annihilation reaction takes place and the simplest one is the steady state.

The first application of diffusion-controlled annihilation to financial markets was suggested by Bak, Paczuski and Shubik to describe market evolution. These authors introduced a series of models based on diffusion-controlled annihilation as a route towards recovering the observed Levy-Pareto “fat-tail” distributions which are said to describe the medium term evolution of the stock market. In BPS model market microstructure, i.e. the structure of reaction zone was not considred at all. For discussions about similarities and differences between BPS and our approaches see . Let us also note that in this paper we do not give references on economical and financial papers on market microstructure. The list of more than half a hundred important publications is given in .

In our model, we shall be studying market microstructure measured over short time scales, so that sociological interactions between traders may be ignored. In the next few sections, we will describe the details of the models (§3.1 and §3.2).

Much research has been done in steady-state diffusion-driven annihilation reactions, especially in physics and chemistry. In applications in physics and chemistry, such reactions are usually in three dimensions (e.g. in a solid or a liquid), but in our model of financial markets, we deal with one dimension only. In the initial paper by Gálfi and Rácz, the properties of the reaction front in a system with segregated initial conditions were studied in the mean field approximation. The mean field approximation, which ignores higher, non-Gaussian corrections, known as fluctuations, works well for dynamics in three dimension, but breaks down in lower dimensions due to the important role of microscopic density fluctuations in one and two dimensions. Numerical simulations were performed in . Analytical calculations by Cardy et. al.  confirmed these numerical results.

This paper is organised as follows. In the next section brief introduction into diffusion driven annihilation reactions (sometimes called DCR, i.e. diffusion controlled reactions) is given. Then applications of DCRs to market microstructure are discussed (mostly based on ) in the section $`3`$. In the last two sections the results of numerical simulations for Minimal and Two-Liquid models are presented. Some of this results confirm earlier known analytical or numerical results, but most of them are new - for example results about bid-offer spread, time to midmarket sale in Minimal model and critical ratio of market order traders in the Two-Liquid model. In the conclusion we discuss obtained results and new interetsing problems worth studying in a future.

## 2 Diffusion driven annihilation reactions

### 2.1 Basic principles of diffusion controlled reactions

Diffusion controlled reactions, or DCRs, are reactions which include as one of the stages the transport of components which can be described by diffusion equations (see for example ). In DCRs, the rate of transport of reacting particles to within reacting distance of one another is the rate-determining step, and the rate of actual chemical combination is assumed to be much faster than that of transport by diffusion. This is a common situation that occurs not only in chemical reactions between atoms, molecules and ions, but also in cases involving electrons, lattice defects, quasi-particles, dislocations, etc.

In classical chemical kinetics, the rate of chemical combination is assumed to be sufficiently low compared with the rate of transport, that the uniform particle distribution is not disturbed and the system of two species may be considered to be in quasi-equilibrium. The rate of reaction is therefore a function of the state of the system, determined by the temperature, pressure and particle concentration. Diffusion space relaxation is fast enough so that diffusion is not the rate-determining step. In the classical approximation, therefore, we may characterize the dynamics by simple functions of state.

This is not the case in the diffusion-controlled regime. Inhomogeneities arise because of local depletion of reactants in regions of high reaction rate. The overall rate of reaction is determined not only by the mean reactant concentration and other functions of state but also by the relative positions and velocities of the particles. The inhomogeneities have linear dimensions of the order of the inter-particle separation. In such cases, the notion of a local concentration is meaningless and the system must be described by a multi-particle distribution function. The search for the probability of particle collisions is much complicated by this. The mutual correlation in the positions of the particles is determined by the chemical reaction, making the probability of collision between two particles a complex function of the coordinator of all the particles in the system. In short, we have come up against the many body problem. Methods of quantum mechanics and statistical mechanics are widely applied in the search for a solution.

DCR theory, which seeks a quantitative understanding of DCRs, consists of two stages. The first stage involves simple models based on single-particle distribution functions. The second stage incorporates the role of correlation effects. There are four basic assumptions:

1. A System containing chemically active ‘primary’ reagents is in a state of thermal equilibrium. In a process consisting of a fast and a slow stage, the fast stage produces species which first achieve a thermal equilibrium, before these unstable species react further by the slow process to reach a chemical, and therefore thermodynamic, equilibrium. This is well illustrated by high energy radiation chemical processes, whereby radiation incident on a sample causes excitations to be produced in the molecules, forming ions, radicals and other reactive species. These species are non-uniformly distributed, usually in the form of tracks or cords. The system first comes to thermal equilibrium through exchange of heat energy between parts at different local temperatures. There then follows a progression towards chemical equilibrium, and thus thermodynamic equilibrium, through the ensuing chemical transformations of the unstable species.
2. Initially, chemically active particles are assumed to be distributed inhomogeneously throughout the volume in the general case.
3. The main postulate of DCR theory is that active particles diffuse according to Fick’s law, $`J=Dn(𝐫)`$ where $`J`$ is the diffusion current, $`D`$ is the diffusion coefficient and $`n(𝐫)`$ is the local active particle concentration. This is applicable in the presence of strong macroscopic concentration gradients. However, in a spatially inhomogeneous distribution, the notion of concentration is meaningless, as the characteristic length scales involved are of the order of the intermolecular separation. Here, we have to resort to analyzing diffusive mobility on a microscopic molecular level. We can imagine the particles situated on an imaginary lattice and hopping from site to site within a short time interval. After each jump, the particle comes into thermal equilibrium with the medium, from assumption 1. The jumps are statistically independent, represented by a random walk. This leads to the use of the theory of stochastic processes.
4. Chemical reactions determine boundary conditions of the differential equations. The system is initially in a state of quasi-equilibrium. The chemical reactions do not change the inhomogeneous state of the system. The problem of how to take into account the limiting transport (diffusive) stages of the chemical transformation was solved by Smoluchowski by considering the probability for a certain particle to remain unreacted by a certain time. The dynamics of unreacted particles are still determined by the diffusion equations. Thus, the probability of reaction can be determined by the flux of mobile particles through the boundary assigned by the selected particle surface. Taking into account chemical transformations in DCR kinetics is reduced to the problem of applying boundary conditions to the diffusion equations.

### 2.2 Applications of diffusion driven annihilation reactions

Diffusion driven annihilation reactions have been well studied in physics and chemistry, partly because of their wide applications across many fields of these subjects. The applications in chemistry are obvious; after all, chemistry is about the study of chemical reactions, a large number of which are diffusion-limited. There are also applications in physics, especially in the field of condensed matter where large numbers of particles are being considered. A newer, and perhaps less obvious, application is in economics, where the market microstructure determined by the short term behaviour of traders may be viewed as a form of diffusion through price space. We shall not go into it in any more detail here as the entire §3 is devoted to this application of DCRs. Here are three examples of applications in chemistry and physics:

1. Atomic species $`A`$ and $`B`$ diffuse through a given volume with diffusion coefficients $`D_A`$ and $`D_B`$. When $`A`$ and $`B`$ are a distance $`R`$ from each other, they recombine in time $`dt`$ with *a priori* probability $`\frac{1}{\tau (R)}dt`$. Thus

$`A+B`$ $``$ $`C`$
$`(\mathrm{inert})`$
in a unidirectional reaction. If we have a bidirectional reaction where the product is not inert,

$$A+BC$$

the reverse process has a decomposition rate of $`\frac{1}{\tau (R)}e^{\beta U}`$ where $`\beta =1/kT`$. Normally, $`U0`$ but if the reaction is mainly that of decomposition, then $`U<0`$.
2. Fluorescence of certain liquid or amorphous semiconductors. The electrons and holes within the semiconductor diffuse separately and independently. When they are in close proximity to one another, they have a finite probability to recombine and emit gamma rays. Here, the unidirectional model applies. It also applies to the diffusion and decay of excitons at recombination (‘scavenger’) sites, an alternative and important mechanism for delayed illumination which involves only a single species.
3. The structure of imperfect solids contains an excess of vacancies (missing atoms) such as those created by radiation damage. Vacancies diffuse throughout the solid until they recombine with interstitial atoms, or diffuse to the surface and effectively disappear. Alternatively, vacancies or interstitials created at the surface can diffuse to the interior. Diffusion-limited aggregation, whereby certain solids grow from vapour or liquid, has proved a popular subject of research.

### 2.3 Uses of quantum field theory in diffusion-limited reactions

Over the past two decades or so, it has been found possible, by many researchers, to use quantum field theory,  to study the motion, diffusion, recombination, and other dynamics of many body systems of *non* quantum mechanical objects, i.e. objects for which $`\mathrm{\Delta }p\mathrm{\Delta }x\mathrm{}`$ plays an insignificant role. It is used as a counting device, and is familiar to many physicists, especially particle and theoretical physicists. In quantum field theory, particles are constrained to execute random walks on the vertices of a space lattice in $`d`$ dimensions. Quantum field theory allows a formal solution of the ‘master equation’ that governs many-body probabilities. To illustrate the method, we describe its application with an example.

The gist of the method may be demonstrated by considering the simplest example of all, radioactive decay. This is a zero dimensional problem, so we can proceed without worrying about diffusion operators. In a group of $`N`$ radioactive atoms, we expect the number remaining undecayed at time $`t`$ to be $`n(t)=N\mathrm{exp}(t/\tau )`$ where $`\tau `$ is the mean lifetime. $`n(t)`$ is the expected number remaining undecayed, averaged over a large ensemble of such $`N`$ atoms. It is known that initial conditions\[17, p981\] affect the short and medium term behaviour of the system, but the asymptotic behaviour for $`t\tau `$ is independent of initial conditions (as borne out by simulations in ).

If $`1/\tau `$ is the rate at which a single nucleus decays, then, for the set of $`n`$ undecayed nuclei, the rate at which a single decay occurs is $`n/\tau `$. The master equation for the probabilities is a linear differential-difference equation,

$$dP(n|t)=\{(n+1)P(n+1|t)nP(n|t)\}\frac{dt}{\tau },$$
(2.1)
which incorporates the two unidirectional processes: decay *into* the state of occupancy $`n`$ from $`n+1`$, and decay *out of* it into $`n1`$. This master equation may be solved by indirection, by associating the state of $`n`$ particles with the $`n`$th excited state of a harmonic oscillator. Consider the harmonic-oscillator raising operator $`a^{}`$ and its conjugate operator $`a`$. The commutation relation $`[a,a^{}]=1`$ holds. If we define $`|0`$ as the ground state with zero occupation number, such that $`a|00`$, the following relations may be obtained. We define $`|n)`$ by

$`a^{}|n)=|n+1),`$ $`a|n)=n|n1).`$
This implies

$`|n)`$ $`=`$ $`(a^{})^n|0,`$
$`a(a^{})^n|0`$ $`=`$ $`n(a^{})^{n1}|0,`$
$`a^{}a(a^{})^n|0`$ $`=`$ $`n(a^{})^n|0.`$
It follows that $`a^{}a`$ is the number operator. We now introduce the crucial concept of a *probability state vector* in which to embed the $`P(n|t)`$:

$$|\mathrm{\Psi }(t))\underset{n=0}{\overset{\mathrm{}}{}}P(n|t)|n)=\underset{n=0}{\overset{\mathrm{}}{}}P(n|t)(a^{})^n|0.$$
(2.2)
The initial condition is that there are $`N`$ atoms at $`t=0`$, so $`P(n|0)=\delta _{n,N}`$, which is equivalent to the probability state $`|\mathrm{\Psi }(0))=(a^{})^N|0`$. The right-hand basis states in this vector space (Fock space) are the $`(a^{})^n|0`$. Given

$$1/n!0|(a)^n(a^{})^n^{}|0=\delta _{n,n^{}}$$

they form a complete orthogonal set. However, the states are not normalized in the conventional manner. Instead, they are normalized with respect to a “reference” state $`(𝒮|0|e^a`$ which is merely a special case of a *Glauber state* $`\alpha 𝒮|0|e^{\alpha a}`$ (itself an eigenvector of the operator $`a^{}`$ with eigenvalue $`\alpha `$).

The norm of any right-hand state $`|\mathrm{\Phi })`$ is defined in terms of $`(𝒮|`$ by the inner product $`(𝒮|\mathrm{\Phi })`$, so that a state is normalized if it satisfies $`(𝒮|\mathrm{\Phi })=1`$. Given that $`(𝒮|a_{}^{}{}_{}{}^{n}=(𝒮|`$ for all $`n0`$, and $`(𝒮|0=1`$, each of the right-hand basis states $`(a^{})^n|0`$ is normalized in this fashion, i.e. $`(𝒮|(a^{})^n|0=1`$ for any $`n0`$. Note that the probabilities remain normalized at all times $`t`$ if they were normalized initially, provided the equations of motion conserve probability (which they always do):

$$(𝒮|\mathrm{\Psi }(t))\underset{n}{}P(n|t)(𝒮|(a^{})^n|0=\underset{n}{}P(n|t)=1.$$

Expectation values may be computed as follows. To evaluate the average number remaining at any time $`t>0`$, we make use of the number operator $`a^{}a`$,

$$(𝒮|a^{}a|\mathrm{\Psi }(t))\underset{n}{}nP(n|t)n.$$

In general, the expectation value of a function of the number operator $`F(a^{}a)`$ may be computed thus:

$$(𝒮|F(a^{}a)|\mathrm{\Psi }(t))\underset{n}{}F(n)P(n|t)F(t).$$

The state vector $`|\mathrm{\Psi }(t))`$ has all the properties of a generating function. First, it satisfies an elementary differential equation equivalent to Eq. (2.1). Second, it yields individual $`P(n|t)`$ by projection onto $`1/n!0|a^n`$. Third, it can be used to obtain the various moments, through contractions with $`(𝒮|`$. The differential equation satisfied by $`|\mathrm{\Psi }(t))`$ is

$$_t|\mathrm{\Psi }(t))=\frac{1}{\tau }\mathrm{\Omega }|\mathrm{\Psi }(t)),$$
(2.3)
in which $`\mathrm{\Omega }=a^{}aa`$, the dimensionless rate operator for this process, is sometimes referred to as the “quantum Hamiltonian” of the model. We illustrate how Eq. (2.3) may be derived from the state vector equation (2.2). The technique is to differentiate Eq. (2.2) with respect to time, substitute in the master equation (2.1) and then try to arrange it into the form of a Schrödinger equation, as in (2.3):

$`_t|\mathrm{\Psi }(t))`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{dP(n|t)}{dt}}|n)`$ (2.4)
$`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\{(n+1)P(n+1|t)nP(n|t)\}{\displaystyle \frac{1}{\tau }}|n)`$
Now, we observe that

$`(a^{}aa)|\mathrm{\Psi })`$ $`=`$ $`(a^{}aa){\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}P(n|t)|n)`$ (2.5)
$`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}nP(n|t)|n){\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}nP(n|t)|n1)`$
$`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}nP(n|t)|n){\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}(n+1)P(n+1|t)|n)`$
which we can identify with the right hand side of Eq. (2.4). Therefore, we can transform the master equation involving probabilities to a Schrödinger-like equation involving $`|\mathrm{\Psi })`$:

$`_t|\mathrm{\Psi }(t))`$ $`=`$ $`{\displaystyle \frac{(a^{}aa)}{\tau }}|\mathrm{\Psi }(t))`$
$`=`$ $`{\displaystyle \frac{1}{\tau }}\mathrm{\Omega }|\mathrm{\Psi }(t))`$
giving us Eq. (2.3). The solution is simply

$$|\mathrm{\Psi }(t))=e^{t/\tau \mathrm{\Omega }}|\mathrm{\Psi }(0))=e^{t/\tau \mathrm{\Omega }}(a^{})^N|0,$$

where the right hand side represents the initial condition of having precisely $`N`$ undecayed nuclei. The exponential in the solution may be expanded in series form and the whole expression substituted back into Eq. (2.3) to check that it satisfies the differential equation.

The above illustration of the method for 0 dimensions can be extended to 1, 2 and 3 dimensions, with the introduction of diffusion and annihilation operators. Although the example shown above is trivial, quantum field theory is a powerful method for solving for the dynamics of particles diffusing in a lattice, and allows analytic solutions to be obtained. It is beyond the scope of this report to discuss further examples of such applications, though the interested reader is referred to ,  for a detailed review.

## 3 The market microstructure of the interdealer broker markets as a reaction zone of one-dimensional diffusion driven annihilation reaction

### 3.1 The minimal model

Because we are only interested in scaling properties, and not in modelling very detailed behaviour, the model is limited to only a few necessary parameters. This allows us to focus attention on the main model feature under examination. In , it was proposed that the feature dominating double auction market microstructure was the great press of numbers of traders, so that the behaviour of, e.g. moments of measures of liquidity, were largely determined by the statistics of the traders. The required model is one in which a large number of individual traders interact according to any plausible dynamcics for the traders, and it is natural to choose the simplest one. Trader dynamics *require* that when a buyer and seller arrive at the same point in price space, they have agreed on a price, and must do a deal and vanish from the price space. And the simplest plausible model of trader motion is a random walk. Since we only wish to model markets in a steady state, we give these traders no net drift.

This model of microstructure of the interdealer broker market can be mapped onto a statistical field theory, in which particles diffuse and annihilate. Thus, we associate to each buyer his bid price, and each seller his offer price, and imagine the two types of trade prices moving around in one-dimensional price space (Fig. 3.1).

This price space is a discrete lattice, with a finite size. The traders, which may be thought of as particles, hop from one lattice site to the next with a certain probability (or they may choose to stay where they are). In this model, each trader trades in a standard size. When a buyer and a seller meet up at the same point in price space, they annihilate, leaving the market (Figs. 3.2 and 3.3).

The reaction may therefore be represented by $`B+S0`$, which is a *diffusion controlled (or diffusion driven) annihilation reaction*. This type of reaction has been studied extensively in the physics and chemistry literature.

We treat the case of a quiescent market, one which is an *approximately steady state* market, and we define our evolving state measure, at any time $`t`$, to be that probability measure which is the steady state solution of our diffusion-annihilation dynamics. As a first approximation, we allow traders to enter only at the ends, as though the buyers start by bidding very low and the sellers offering very high. One can solve this model by numerical simulation. However, because of its similarity to many well-known analytically soluble models (for example, in ), we can go much further. Approximate analytical methods are available from Cardy et al.

### 3.2 The two-liquid model

The minimal model is a suitable approximation of the quiescent periods of a market’s behaviour, but when it is in the process of a long crash, it does not look quiescent and the model breaks down. It is not sufficient merely to cause an imbalance in the rates of insertion of buyers and sellers, because such a method would not reproduce the well-known phenomenon of the widening of the bid-offer spread. It is believed that the missing element is market order traders. During a crash, most traders, desperate to sell, attempt to hit the best bid directly, as a market order. In such a scenario, although limit order trading continues as before, the large proportion of market order trades cannot be neglected and indeed they have a significant effect on market dynamics. In the two-liquid model, we introduce a second type of trader, the market order trader, who tries to hit the best bid/offer directly, but without displaying his/her bid/offer price on the trading screen. These market order traders trade with limit order traders, who display their prices on the trading screen, but not with one another. If limit order traders are represented by $`B`$ and $`S`$ and market order traders by $`B^{}`$ and $`S^{}`$, then the reactions $`B+S0`$, $`B^{}+S0`$, and $`B+S^{}0`$ can occur, but *not* $`B^{}+S^{}0`$, since market order traders are invisible to one another. This model is able to reproduce the familiar widening of the bid-offer spread during a crash.

### 3.3 Dimensional analysis

Because we have relatively few degrees of freedom (parameters) in our model, we can use dimensional analysis to obtain the scaling forms of many market quantities, as a check on intuition. Starting from the operator evolution equation, we may differentiate it to yield an operator differential equation. This partial differential equation consists of one coefficient only, namely the diffusion coefficient $`D`$, with dimensions $`x^2/t`$, i.e. $`D(\mathrm{dollar})^2/\mathrm{sec}`$ (see §A.1 for derivation of the diffusion coefficient). In addition, there is another dimensionful quantity in the boundary conditions, $`J`$, the rate of insertion of traders at the boundaries, which is equal to the deal rate in a steady state market, with dimensions $`J1/(\mathrm{sec})`$. From $`D`$ and $`J`$ we can construct quantities with any dimensions containing (dollar) and (sec), since $`(\mathrm{dollar})\sqrt{D/J}`$ and $`(\mathrm{sec})1/J`$. Therefore, the expectation value $`X`$ of a general quantity $`X`$ with dimensions $`[X]=(\mathrm{dollar})^\mathrm{m}/\mathrm{sec}^\mathrm{n}`$ must be proportional to $`(D/J)^{m/2}J^nD^{m/2}J^{nm/2}`$.

So far, we have ignored, in our dimensional analysis, of the lengths scales $`L`$, the size of price space, and $`\delta S`$, the lattice spacing. These two lengths represent the upper and lower length scales in our model. We are justified in doing so by the fact that if we imagine a Fourier analysis of the dynamics, most of the ‘action’ would be taking place away from these extreme length scales, in the intermediate region between the two. Ultimately, the justification comes from numerical simulations. However, these length scales do add logarithmic corrections to the scaling laws, as shown by Cardy et al. This is described in §3.6.

### 3.4 Market parameters

The model allows us to calculate the relationships described in the introduction, but also many more. In general, whatever initial conditions we start with, in the asymptotic limit, we will reach a steady-state, with a reaction front, where the buyers meet the sellers. At the centre of the reaction front is the best bid, best offer and midmarket. Beyond the best bid/best offer, we expect to see an increasing density of buyers/sellers. Traders are inserted at the edges at a rate equal to the deal rate.

The statistics describing the state of the system can be read off a trading screen. The best bid locates the top of the lower edge and the best offer the bottom of the upper edge (see Figs. 3.1 and 3.3). Although there is a large number of interesting parameters, in this project we shall be restricting our attention to a small subset of them: best bid $`B(t)`$, best offer $`O(t)`$, bid-offer spread $`Spr(t)O(t)B(t)`$, and midmarket $`M(t)(1/2)(B(t)+O(t))`$. Note that $`B(t)`$, $`O(t)`$, $`Spr(t)`$ and $`M(t)`$ are not Markov random variables, and do not satisfy a stochastic differential equation, because they are subject to jumps. In addition to the above parameters, we shall also be interested in recording the instantaneous deal rate $`DR`$ and the total number of traders in the market $`NUM`$.

### 3.5 The scaling laws

Financial markets, like many other systems, exhibit scaling laws which are clearly observed by practitioners in the field. These scaling laws are related to the concept of liquidity, which is measured in several different ways. It may be measured by the deal rate $`J`$, the time to midmarket sale, the bid-offer spread, and the trader densities near the best bid/offer. Here we shall concentrate on several scaling laws which were introduced in .

* Bid Offer Spread

$$Spr=O(t)B(t)$$
(3.1)
* Trade price variance (or midmarket variance)

$$w^2=M^2(t)M(t)^2$$
(3.2)
* Fluctuations in the Bid-Offer Spread

$$(\mathrm{\Delta }Spr)^2=\mathrm{var}(Spr)=(O(t)B(t))^2O(t)B(t)^2$$
(3.3)
* Instantaneous Dealing Time, and its Fluctuations

We can define $`\tau _S`$ as the time between successive trades. If after a time, more than one trade occurs in the same time step, then, in order to have a valid continuum limit for the quantity, we must consider the time between the multiple trades themselves to be zero. We may visualize it thus, assuming time is continuous: one trade occurs at $`t=0`$, and then no trades occur for some time. The next trade occurs at $`t=9`$, and the next at $`t=9+\delta `$. These are distinct times, if we use continuous time, and we would say that the two values of $`\tau _S`$ are 9 and $`\delta `$. If we now allow $`\delta `$ to tend to zero, we have 9 and 0 as the values of $`\tau _S`$, but we see that such a configuration is equivalent to the discrete model of time, where we have no trade for 9 time steps followed by 2 at the same time. Thus, the correct way to deal with “simultaneous” annihilations is to treat them as separate ones, separated by zero time.

We can define its fluctuation as

$$\mathrm{\Delta }\tau _{S}^{}{}_{}{}^{2}=\mathrm{var}(\tau _S)=\tau _{S}^{}{}_{}{}^{2}\tau _S^2$$
(3.4)
* Asymmetric fluxes

It is possible to consider the minimal model with different fluxes for the buyers and the sellers. In  the following theoretical result, valid for “small” times, was found:

$$x_0(t)\frac{2D\mathrm{\Delta }J}{\overline{J}L}t,$$
(3.5)
where $`x_0(t)`$ is the position of the midmarket, which varies with time $`t`$. The flux of sellers from the upper boundary is $`\overline{J}+\mathrm{\Delta }J`$ whilst the corresponding flux of buyers from the lower boundary is $`\overline{J}\mathrm{\Delta }J`$.

We expect $`Spr(J)`$ and $`w^2`$ to tend to 0 as $`J\mathrm{}`$. From dimensional analysis, $`Spr`$ and $`w`$ have dimensions of *dollars*, so $`Spr(J)\sqrt{(}D/J)`$, $`w(J)\sqrt{(}D/J)`$ consistent with intuition. Using the results of , we see that there is a logarithmic correction for $`w`$. It is also possible to measure fluctuations in the spread, $`\mathrm{\Delta }Spr`$ or $`\mathrm{var}(Spr)`$. This quantity should have the same scaling law as $`Spr`$, since it has the same dimensions. Time to midmarket sale $`\tau _S`$ is expected to go as $`1/J`$, from dimensional analysis; in fact, its average should be equal to $`1/J`$. Similarly, its fluctuation, $`\mathrm{\Delta }\tau _S`$ should have the same scaling law. There are many other scaling laws which we will not go into here as they will not be tested later; for more information, see .

### 3.6 Analytical results

Although numerical results are more precise, it is desirable to work out some scaling laws for the model using analytical means, as these methods supplement numerical simulations by providing intuition. Cardy et. al. has developed an approximation scheme known as mean field theory, replacing the evolution equations for the operators with a partial differential equation for the density configuration with the greatest probability mass, ignoring the “fluctuations”. There is no system of higher order corrections or error estimate with this method. In spite of its limitations, its predictions have been found to coincide with numerical simulations . It is beyond the scope of this paper to describe in detail the derivation of the analytical results, so we shall simply state the most important result for our simulation—the logarithmic correction to the scaling law for $`w`$:

$$w=\left[\frac{\mathrm{ln}(cL/w)}{\pi (J/D)}\right]^{1/2}$$

This was derived by the addition of an effective noise term to the differential equation satisfied by the density difference of buyers and sellers. However the numerical value of constant $`c`$ is unknown from analytical calculations and we shall find it later from numerical data.

The second analytic result we shall be using relates to the speed of the moving midmarket in the case of asymmetric fluxes of buyers and sellers. The result here is (for more details see page 34 in )

$$x_0(t)\frac{2D\mathrm{\Delta }J}{\overline{J}L}t$$

and was derived using a time-dependent non-stochastic part in the density difference equation. It is only valid for small times, such that $`T<L^2/D`$.

## 4 Computer simulations of the minimal model

### 4.1 Description of simulation

We shall start from a Monte-Carlo simulation of the minimal model.

The simulation starts by randomly inserting traders into either half of price space to set up a random initial configuration. The main loop is entered. Traders are inserted at the edges, and then each type of trader (buyer and seller) is allowed to diffuse separately and annihilate. Each trader hops to the left or the right at every time step with probability 0.5 ($`p=1`$). This leads to a diffusion coefficient of 1/2, because $`D=a^2/2\tau `$ (see §A.1), and $`a=\tau =1`$ in our units. The annihilation routine is called after each diffusion routine to ensure that pairs of mutually annihilating traders do not cross over each other. This is in accordance with the procedure outlined in the appendix of , where the time-evolution operator is constructed from diffusion and annihilation operators, with annihilation occurring after each diffusion of buyers and sellers (cf.  which has a slightly different way of doing it). The tasks performed within each time quantum $`\tau `$ may be summarized as:

1. Diffuse buyers
2. Annihilate overlapping traders
3. Diffuse sellers
4. Annihilate overlapping traders
5. Call `getstats()` routine to store market parameters

This method we shall call MINIMAL. At the end of each time step, the `getstats()` routine is called, which obtains interesting statistics from the state variable, such as best bid, best offer, bid-offer spread, midmarket, instantaneous deal rate, etc. Thus the diffusion-annihilation process is repeated over and over again. Now, not all the data generated at each time step is recorded, simply because it would be impractical and unwieldy to have a set of two million results to deal with. Instead, results such as bid-offer spreads and midmarkets are averaged over a large number `DISP_INT` of time steps. In the end, we would like to take the average<sup>§</sup><sup>§</sup>§Note that the average taken of the results generated by computer simulation is a *time* average, whilst the “average” used in the scaling laws and the analytical results ($`\mathrm{}`$) is an *ensemble* average. These two are generally not the same. However, for a stationary ensemble (i.e. one that has no preferred origin in time), the ergodic assumption, that the system will in the course of a sufficiently long time pass through all the states accessible to it, leads to the fact that the time average is equal to the ensemble average. The financial market, when in equilibrium, is well-approximated by a stationary ensemble, and so the equality of the ensemble and time averages applies. That is why it is acceptable to compare theoretical ensemble averages with the time averages obtained from numerical simulations, which are the same as ensemble averages. See  for proof. and variance of all the results, not just the results averaged over `DISP_INT` steps. At first, one might be inclined to think that the process of averaging over some steps causes one to lose information. In fact, it is possible, by simultaneously keeping track of the mean of the squares of the quantities, to allow the mean and variance of the entire ensemble to be recovered, as though we had kept track of every single result (for proof see Appendix §A.2).

The data produced from the program were recorded in a data file. Note that in the simulation, the quantities considered were dimensionless, and may be defined as follows:

$`T=N_1\tau `$
$`L=N_2a`$
$`J\tau =N_3`$
where $`\tau `$ is the time between each step in the simulation (the time quantum), $`T`$ is the total time of the simulation, $`L`$ is the width of price space, $`a`$ is the lattice spacing of price space, and $`J`$ is the rate of insertion of traders (no. of traders per second). $`N_1`$, $`N_2`$, and $`N_3`$ may be thought of as the dimensionless counterparts of $`T`$, $`L`$ and $`J`$, respectively. We may simply choose our length and time units such that $`a=\tau =1`$, so that lengths are measured in units of the length quantum $`a`$ and time intervals in units of the time quantum $`\tau `$. In these units, $`N_1=T`$, $`N_2=L`$ and $`N_3=J`$, and so we may revert to the natural units, which have been made dimensionless by a judicious choice of units.

### 4.2 Preliminary results for $`L=200`$ using MINIMAL

The simulation of the minimal model was initially performed for 100,000 time steps, with averaged results recorded every 100 or 1000 steps. The number of steps in the simulation was limited by the length of time taken to complete it: for each run of 100,000 steps and one particular value of $`J`$, the run-time would be between 1 and 2 hours. The time taken to obtain an entire set of results was a few days. There were 200 points in price space ($`L=200`$), and the initial condition consisted of 100 buyers and 100 sellers distributed randomly on either side of price space, with a gap of 20 empty price points in the middle. ¿From this condition, the system was allowed to evolve according to the diffusion and annihilation laws outlined earlier. Every `DISP_INT` time steps, the following averaged quantities were recorded: best bid ($`B`$), best offer ($`O`$), best bid size ($`BS`$), best offer size ($`OS`$), midmarket ($`M`$), bid-offer spread ($`Spr`$), deal rate ($`DR`$), and total number of traders in the market-place ($`NUM`$).

It was found that the system came to equilibrium after about 50,000 time steps. The total number of traders in the market ($`NUM`$) was used as a suitable measure of the state of equilibrium of the system, because it was discovered that $`NUM`$ grew with time and tended asymptotically to an equilibrium value (Fig. 4.1). We have also included a picture of how the bid, offer and midmarket vary with time in the simulation (Fig. 4.2).

Thus, averages of quantities were only taken over the range of time steps for which the system was in equilibrium, in this case, the last 50,000.

The scaling laws we are interested in are the following:

$`w^2`$ $`=`$ $`{\displaystyle \frac{D}{\pi J}}\mathrm{ln}\left({\displaystyle \frac{cL}{w}}\right),`$ (4.1)
$`Spr`$ $``$ $`\sqrt{D/J},`$ (4.2)
$`\mathrm{var}(Spr)`$ $``$ $`\sqrt{D/J}.`$ (4.3)
To this end, we plot graphs of $`\mathrm{ln}w^2`$ against $`\mathrm{ln}J`$, and $`\mathrm{ln}Spr`$ against $`\mathrm{ln}J`$. We can see that each run of the simulation produces *one* point on a graph, since we average over a long time to obtain mean values for our market parameters. To obtain any reasonable graph requires at least 10 points or so, hence the long times taken to obtain results.

The results obtained (Figures 4.6, 4.6, 4.6 and 4.6) were of a poor quality and did not yield the expected straight lines. The values of $`J`$ used were 1, 2, 4, 8, 16, 32 and 64.

One can see that these values of $`J`$ are too high and that in this regime, the market is very dense and effectively frozen. This is not the regime described by the minimal model, and it is not surprising that the results do not agree with the predictions. From the data, one can see that $`w<a`$, $`\mathrm{var}(Spr)<a`$, and $`Spra`$. The reason for the poor results is the following: the scaling laws, derived using dimensional analysis, relied on the fact that the length scales of meaningful quantities were far from $`L`$ and $`a`$. Thus, for dimensional analysis to hold, and the scaling laws obtained thereby to be valid, we require that

$$aw,Spr,\mathrm{var}(Spr)L.$$
(4.4)
This means that $`J`$ must be much less than 1, since all the length quantities are of order $`\sqrt{(}D/J)`$.

### 4.3 Successive improvements to preliminary results

We outline a sequence of improvements to the simulations as follows:

| Method | $`L`$ | $`T`$ | Range of $`J`$ |
| --- | --- | --- | --- |
| MINIMAL | 200 | $`1\times 10^5`$ | $`J>1`$ |
| MINIMAL | 200 | $`1\times 10^5`$ | $`J<1`$ |
| MINIMAL | 200 | $`1\times 10^6`$ | $`J<1`$ |
| MINIMAL | 1000 | $`2\times 10^6`$ | $`J<1`$ |

The average time taken to obtain each set of results was about 3 or 4 days. The simulations were repeated for fractional $`J`$, from 0.001 to 0.1, for $`L=200`$. The results obtained (MINIMAL, $`L=200`$, $`T=1\times 10^5`$) were better than before, and the points lay more closely on a straight line. As the simulation time is proportional to $`J`$ (see §4.13), it becomes possible, with $`J`$ reduced by several orders of magnitude, to increase the number of time steps per run to 1 million whilst keeping the total duration reasonable, with the statistics being taken over the last 950,000 steps. Thus, it was decided to do the simulations again, but with 1 million steps. The results obtained (MINIMAL, $`L=200`$, $`T=1\times 10^6`$) were even better, though the graphs were still not very convincing, especially the one of $`Jw^2`$ against $`\mathrm{ln}w`$, which should have been a straight line, according to Eq. (4.1), but in fact the points were so scattered that it was difficult to draw any definite conclusion.

The smallness of $`L`$ (200) restricted the choice of $`J`$ too much: it was desired to keep $`w`$ at least a factor of 5 away from either $`a`$ or $`L`$, so the range of allowed $`w`$ was from 5 to 40, which was just one order of magnitude. Thus, it was resolved to increase $`L`$ to 1000, which expanded the range of permitted $`w`$ to $`5<w<200`$. However, the increase in the size of the system caused a corresponding increase in the time taken for the system (especially trader number $`NUM`$) to reach equilibrium. It was noticed that quasi-equilibrium was not attained until after about 1 million steps. The duration of the simulation was therefore increased to 2 million steps, with only the last 1 million steps used for obtaining statistics.

### 4.4 Encouraging results using MINIMAL for $`L=1000`$ and $`T=2\times 10^6`$

Many different values of $`J`$ were used, all of which produced values of $`w`$ and $`Spr`$ roughly within the range specified above, with most interesting lengths being at least a factor of 5 away from either $`a`$ or $`L`$. Here, the initial condition consisted of 50 buyers and 50 sellers, separated by a gap of 20. One of the things learned from these simulations (obtainable by intuition, and supported by ) is that the asymptotic behaviour of the market does not depend on the initial conditions.

38 different values of $`J`$ were used, producing a set of 38 points on the graphs (Figs. 4.7, 4.8, 4.9 and 4.10).

It was noticed that for the first three graphs at least, the points fit a line better towards the bottom right of the graph, while the points towards the top left tended to be more scattered. It was likely that the points towards the top left (small $`J`$ and large $`w^2`$) were less accurate because of the possibility of the reaction front crashing into either edge of price space. In addition, for graphs of $`\mathrm{ln}w^2`$, $`\mathrm{ln}Spr`$ and $`\mathrm{ln}\mathrm{var}(Spr)`$ against $`\mathrm{ln}J`$, the analysis in §4.5 requires that $`\mathrm{ln}w\mathrm{ln}cL`$ be satisfied, and this is not satisfied for the points towards the top left, for which $`\mathrm{ln}w5`$. Therefore, we must choose a small section of the bottom right of the graph and obtain the gradient and intercept from that.

Thus, it was decided to take only the last 20 points and fit a line through them, in the hope of acquiring more accurate statistics (Figs. 4.11, 4.12, 4.13 and 4.14).

Here is a summary of the results, for the first set of data (all 38 points):

| Figure | Abscissa | Ordinate | Gradient | Intercept | $`R^2`$ |
| --- | --- | --- | --- | --- | --- |
| 4.7 | $`\mathrm{ln}J`$ | $`\mathrm{ln}w^2`$ | $`0.814\pm 0.015`$ | $`0.41\pm 0.26`$ | 0.988 |
| 4.8 | $`\mathrm{ln}J`$ | $`\mathrm{ln}Spr^2`$ | $`0.925\pm 0.014`$ | $`0.903\pm 0.24`$ | 0.992 |
| 4.9 | $`\mathrm{ln}J`$ | $`\mathrm{ln}\mathrm{var}(Spr)`$ | $`0.924\pm 0.014`$ | $`0.38\pm 0.25`$ | 0.991 |
| 4.10 | $`\mathrm{ln}w`$ | $`Jw^2`$ | $`0.174\pm 0.013`$ | $`1.006\pm 0.089`$ | 0.844 |

and here is the same for the second set of data (last 20 points):

| Figure | Abscissa | Ordinate | Gradient | Intercept | $`R^2`$ |
| --- | --- | --- | --- | --- | --- |
| 4.11 | $`\mathrm{ln}J`$ | $`\mathrm{ln}w^2`$ | $`0.919\pm 0.016`$ | $`0.052\pm 0.1`$ | 0.995 |
| 4.12 | $`\mathrm{ln}J`$ | $`\mathrm{ln}Spr^2`$ | $`0.9904\pm 0.0024`$ | $`0.592\pm 0.017`$ | 0.99989 |
| 4.13 | $`\mathrm{ln}J`$ | $`\mathrm{ln}\mathrm{var}(Spr)`$ | $`1.001\pm 0.003`$ | $`0.744\pm 0.023`$ | 0.99980 |
| 4.14 | $`\mathrm{ln}w`$ | $`Jw^2`$ | $`0.107\pm 0.024`$ | $`0.888\pm 0.076`$ | 0.523 |

It must be noted that for all the graphs, with the exception of $`Jw^2`$ against $`\mathrm{ln}w`$, the correlation coefficients are all 0.98 or above, which means that the results strongly support a scaling relationship of some sort, the exact details of which we shall now discuss.

### 4.5 The $`w^2`$ scaling law

Now that we have some reasonably detailed results, we may begin our analysis of them in earnest. The fact that the plot of $`\mathrm{ln}w^2`$ against $`\mathrm{ln}J`$ is a straight line with a gradient close to –1 suggests that $`w^2`$ is inversely proportional to $`J`$, which is correct according to Eq. (4.1) if we ignore the logarithmic correction. This is simply a consequence of the dimensional analysis in §3.3.

Next, we seek to observe the logarithmic correction. Consider the plot of $`\mathrm{ln}w^2`$ against $`\mathrm{ln}J`$. From Eq. (4.1), and from $`D=a^2/2\tau =1/2`$, we can write

$`\mathrm{ln}w^2`$ $`=`$ $`\mathrm{ln}J+\mathrm{ln}\left[{\displaystyle \frac{1}{2\pi }}\mathrm{ln}\left({\displaystyle \frac{cL}{w}}\right)\right]`$
$`=`$ $`\mathrm{ln}J+\mathrm{ln}\left\{{\displaystyle \frac{1}{2\pi }}\left[\mathrm{ln}(cL)\mathrm{ln}w\right]\right\}`$
$`=`$ $`\mathrm{ln}J+\mathrm{ln}\left\{{\displaystyle \frac{1}{2\pi }}\mathrm{ln}(cL)\left[1{\displaystyle \frac{\mathrm{ln}w}{\mathrm{ln}cL}}\right]\right\}`$
$`=`$ $`\mathrm{ln}J+\mathrm{ln}\left[{\displaystyle \frac{1}{2\pi }}\mathrm{ln}(cL)\right]+\mathrm{ln}\left[1{\displaystyle \frac{\mathrm{ln}w}{\mathrm{ln}cL}}\right].`$
Here, we make the approximation $`\mathrm{ln}(1x)x`$, ignoring the higher terms $`x^2/2x^3/3\mathrm{}`$, which is valid for $`x1`$, on the third term on the right hand side. In this case, the approximation is valid if $`\mathrm{ln}w\mathrm{ln}cL`$. Letting $`\mathrm{ln}w=1/2\mathrm{ln}w^2`$, and grouping terms of $`\mathrm{ln}w^2`$ on the left hand side, we obtain,

$`\left(1+{\displaystyle \frac{1}{2\mathrm{ln}cL}}\right)\mathrm{ln}w^2`$ $`=`$ $`\mathrm{ln}J+\mathrm{ln}\left({\displaystyle \frac{\mathrm{ln}cL}{2\pi }}\right).`$
For $`2\mathrm{ln}cL1`$, we can make the further approximation of $`\left(1+\frac{1}{2\mathrm{ln}cL}\right)^11\frac{1}{2\mathrm{ln}cL}`$, yielding

$`\mathrm{ln}w^2`$ $`=`$ $`\left(1{\displaystyle \frac{1}{2\mathrm{ln}cL}}\right)\mathrm{ln}J+\left(1{\displaystyle \frac{1}{2\mathrm{ln}cL}}\right)\mathrm{ln}\left({\displaystyle \frac{\mathrm{ln}cL}{2\pi }}\right)`$ (4.5)
$`=`$ $`\alpha \mathrm{ln}J+\alpha \mathrm{ln}\left({\displaystyle \frac{\mathrm{ln}cL}{2\pi }}\right),`$ (4.6)
which defines $`\alpha `$. If we assume that $`cO(1)`$ and take $`L=1000`$, then the logarithm in $`\alpha `$ is positive, making $`\alpha `$ slightly smaller than unity. From measuring $`\alpha `$, we may therefore deduce the value of $`c`$. In this way, it is possible to see the subtle effect of the logarithmic correction. Note that the second approximation is not strictly necessary and one can use the full expression for $`\alpha `$ if one wishes, and this is indeed what we shall do.

We make use of the graph of the last 20 points because we want the approximation $`\mathrm{ln}w\mathrm{ln}cL`$ to hold, as otherwise the foregoing analysis would be invalid. From Fig. 4.11, and using the full expression for $`\alpha `$ (i.e. without making the binomial expansion for the reciprocal) we find that the gradient of the graph is $`0.919\pm 0.016`$, making $`\alpha =0.919`$ and $`c=0.29`$, but because $`c`$ is inside a logarithm, the uncertainty in $`\alpha `$ is magnified in $`c`$, so that it lies in range $`0.11<c<1.33`$.

We can also obtain a value for $`c`$ from the intercept, which is $`0.052\pm 0.1`$, another wildly imprecise quantity. This method of calculating $`c`$ gives 0.38, with limits $`0.20<c<0.75`$. The intercept appears, therefore, to provide a more precise determination of $`c`$, though, as we shall see later, that is not always the case.

There is another technique one can try in extracting information from the graphs. Instead of allowing only $`c`$ to be determined from the graphs, we may attempt to determine the constant $`1/\pi `$ from the gradient and the intercept, since each graph has two degrees of freedom. Therefore, we replace $`1/\pi `$ with a variable, say, $`\mathrm{\Omega }`$, and the equation becomes

$$\mathrm{ln}w^2=\alpha \mathrm{ln}J+\alpha \mathrm{ln}\left(\frac{\mathrm{\Omega }}{2}\mathrm{ln}cL\right),$$

and we may determine $`\mathrm{\Omega }`$ using the data. Here, we have $`\mathrm{\Omega }0.333`$, with the uncertainty $`0.236<\mathrm{\Omega }<0.453`$. Note that, the value predicted from theory is $`1/\pi =0.3183`$, which lies well within the experimental error. This is quite encouraging as it agrees with the result derived theoretically (in spite of its 33% uncertainty!)

It is possible, for the sake of comparison with §4.6 and §4.7, to analyze the data for $`w^2`$ as though its scaling law had no logarithmic correction. It is not really appropriate, because of the relatively pronounced (–0.919) deviation from –1, but we shall quickly do it for subsequent comparison. Assuming that $`w^2`$ is proportional to $`D/J`$ with constant of coefficient $`\lambda `$,

$$w^2=\frac{\lambda D}{J},$$

we can take the log of both sides

$$\mathrm{ln}w^2=\mathrm{ln}\lambda D\mathrm{ln}J$$

and calculate $`\lambda `$ from the intercept, knowing that $`D=1/2`$. The intercept is $`0.052\pm 0.1`$, giving us $`\lambda =1.90`$, within the limits $`1.72<\lambda <2.10`$, or $`1.9\pm 0.2`$.

We have another graph from which we may obtain information: the graph of $`Jw^2`$ against $`\mathrm{ln}w`$. For this graph, it is prudent to take the one for all 38 points, since no approximations have been made in the scaling law we are trying to verify, and the more points we have on the graph, the more confident we may be of its statistics. It is simply

$$Jw^2=\frac{\mathrm{\Omega }}{2}(\mathrm{ln}cL\mathrm{ln}w)$$
(4.7)
which comes from a trivial re-arrangement of Eq. (4.1). We have once again allowed $`1/\pi `$ be determined from the graph. The gradient is $`0.174\pm 0.013`$, from which we may infer $`\mathrm{\Omega }0.348\pm 0.026`$. This is close to the predicted $`1/\pi =0.318`$, although it actually lies outside the range of experimental error (the discrepancy is 9%). We can, with reasonable confidence, now set $`\mathrm{\Omega }`$ to its theoretical value, $`1/\pi `$. From the intercept, $`1.006\pm 0.089`$, we may infer $`c=0.556`$ with the limits being $`0.32<c<0.97`$. These limits are slightly different from those found from the first graph (Fig. 4.11), but not by much. Both methods give the same order of magnitude estimate $`c0.5`$.

### 4.6 The $`Spr`$ scaling law

We shall make use of Figs. 4.8 and 4.12, the latter consisting of the last 20 points of the former. It is interesting to notice that the correlation coefficients for both graphs are very high, and indeed, that of the second graph is so close to 1 that it was rounded up to 1 by the spreadsheet that produced the graph. Once again, the straight lines with gradients so close to –1 confirm the approximate $`Spr\sqrt{(}D/J)`$ law. It is not known (see ) whether there is a logarithmic correction to $`Spr`$ or not. If there is, one might expect it to be of the form

$$Spr^2=\frac{\sigma D}{J}\mathrm{ln}\left(\frac{sL}{Spr}\right),$$

in analogy with Eq. (4.1). We have left the constant in front of $`D/J`$ as a parameter ($`\sigma `$) to be determined.

To look for the logarithmic correction, we can perform a similar analysis as we did in §4.5, making the approximation $`\mathrm{ln}Spr\mathrm{ln}cL`$, to obtain

$$\mathrm{ln}Spr^2=\alpha \mathrm{ln}J+\alpha \mathrm{ln}\left(\frac{\sigma }{2}\mathrm{ln}sL\right),$$

where

$$\frac{1}{\alpha }=1+\frac{1}{2\mathrm{ln}sL}.$$

To ensure the validity of the approximation, we shall use the graph of the last 20 points only (Fig. 4.12). There is another reason for putting more faith in the statistics yielded by the second graph: if one looks carefully along the points of the first graph, one sees that the points towards the lower right of the graph lie along a line that does not quite coincide with the line of best fit constructed from all of the points. This is because of some deviation from linearity towards the top left of the graph. It would be wise, therefore, to take only the results towards the bottom right, as in Fig. 4.12. The gradient is $`0.9904\pm 0.0024`$, which would imply $`s2.5\times 10^{19}`$, $`7.6\times 10^{14}<s<8.8\times 10^{26}`$. The uncertainty covers 12 orders of magnitude! The intercept is $`0.592\pm 0.017`$, giving us $`\sigma 0.07`$, with limits $`0.054<\sigma <0.087`$. We shall try to understand the meaning of the values obtained. Consider a rearrangement of the scaling law

$$Spr^2=\frac{\sigma D}{J}\left[\mathrm{ln}s+\mathrm{ln}\left(\frac{L}{Spr}\right)\right],$$

The largest possible value for $`\mathrm{ln}(L/Spr)`$ occurs when $`Spr`$ is at its smallest, around 5. Thus, the maximum value for the second logarithm is 5.3 (minimum is about 1.6). The first logarithm, however, is $`\mathrm{ln}s`$, which is 44 ($`34<\mathrm{ln}s<62`$). Therefore, the right hand side is dominated completely by the large constant $`\mathrm{ln}s`$, and the logarithm involving $`L`$ has almost no effect on the scaling law. This is a clear sign that there is no logarithmic correction, since the variation of the second logarithm has negligible effect on the variation of the *sum* of the two logarithms, on which the scaling law depends. This can be approximated by ignoring the variation of the second logarithm with $`Spr`$, and taking 3 as the approximate average value for it, leaving $`\mathrm{ln}s47`$, with limits $`37<\mathrm{ln}s<65`$.

$$Spr^2=\frac{\sigma D}{J}\mathrm{ln}s,$$

which is,

$$Spr^2=\frac{\mu D}{J},$$

where the coefficient $`\mu =\sigma \mathrm{ln}s3.3`$.

Now, if we assume that there is no logarithmic correction, we can try to extract from the data the constant of proportionality for the scaling law

$$Spr^2=\frac{\mu D}{J}.$$

Taking logs, we have

$$\mathrm{ln}Spr^2=\mathrm{ln}(\mu D)\mathrm{ln}J.$$

The gradient agrees well with the predicted value of –1. We know that the diffusion coefficient $`D`$ is 1/2, so from the intercept, $`\mu 3.615`$, with limits $`3.55<\mu <3.68`$, or $`3.615\pm 0.065`$. Thus, the two methods, one assuming a logarithmic correction, and the other not, produce results which agree with one another (the 3.3 obtained earlier has an associated uncertainty which is rather difficult to work out, so we will not do that here).

### 4.7 The $`\mathrm{var}(Spr)`$ scaling law

We look at Figs. 4.9 and 4.13. Both have extremely high correlations (almost perfect correlation for Fig. 4.13), so we can have confidence in the results. The results are unambiguously linear, confirming the scaling relationship $`\mathrm{var}(Spr)\sqrt{(}D/J)`$. The only thing remaining is the determination of the existence of logarithmic correction.

The gradient of the graph is $`1.001\pm 0.003`$, which puts the exact value of –1 within the range of experimental error. Even if there is some deviation from –1, it is minute. This is strong evidence *against* the existence of logarithmic correction for $`\mathrm{var}(Spr)`$. There is no point in trying to apply the same analysis as before to find out the constant inside the logarithm, because here, the gradient can lie on either side of –1, causing the value of that constant to be anywhere between the exponential of a large negative number and the exponential of a large positive number, i.e. many orders of magnitude. The fact that the constant can assume these extreme values implies that the logarithm of the constant would dominate the scaling law, when it is either very small, in which case the logarithm of the constant would be a large negative number, or very large, in which case the log would be a large positive number. It is fair to conclude that $`Spr`$ appears from the data *not* to have any logarithmic correction.

We may proceed to find the constant of proportionality, $`\nu `$, defined by

$$\mathrm{var}(Spr)=\frac{\nu D}{J}.$$

Taking logs,

$$\mathrm{ln}\mathrm{var}(Spr)=\mathrm{ln}(\nu D)\mathrm{ln}J.$$

The gradient has already been verified to be –1. The intercept is $`0.744\pm 0.023`$. Therefore, $`\nu 0.95`$, subject to $`0.93<\nu <0.97`$, or $`0.95\pm 0.02`$.

### 4.8 Calculation of $`c`$ for different values of L

¿From Eq. (4.6), we found that plotting graphs of $`\mathrm{ln}w^2`$ against $`\mathrm{ln}J`$ allowed the logarithmic parameter $`c`$ to be determined. After further discussion, it was decided to investigate this parameter and try to see whether and how it varies with $`L`$. To this end, simulations were run for different values of $`L`$, within the correct regime for $`J`$, such that $`aJL`$. The observance of this regime was now especially important, owing to the $`\mathrm{ln}(1x)x`$ approximation we made in our derivation of Eq. (4.6).

The results obtained were rather poor. Correlations were low and the value of $`c`$ did not appear to show any kind of relationship to the value of $`L`$. We shall not present these unhelpful results here. Instead, we shall describe a new method of simulation, which will be used to obtain better results.

### 4.9 MINIMAL1: a new method

If we take a detailed look at the variation of something like $`Spr`$ with time, it will be seen that the MINIMAL method of simulation, which we have been using from the outset, does *not* allow $`Spr`$ to go to zero. This is because of the sequence of steps undertaken at each time quantum $`\tau `$. Only after both species of traders have diffused and annihilated do we take measurements of the system using `getstats()`. In practice, annihilations in the marketplace are not instantaneous, and for some of the time at least, overlapping buyers and sellers *can* exist. This suggests an alternative method of simulation, MINIMAL1, based on the following 8 steps:

1. Diffuse buyers
2. Call `getstats()` routine—overlapping traders possible
3. Annihilate overlapping traders
4. Call `getstats()` routine—no overlapping traders
5. Diffuse sellers
6. Call `getstats()` routine—overlapping traders possible
7. Annihilate overlapping traders
8. Call `getstats()` routine—no overlapping traders

This method represents a more detailed observation of the market as it includes intermediate overlapped states which were ignored in MINIMAL. MINIMAL1 allows the bid-offer spread to vanish just before an annihilation, and can be said to produce a more accurate profile of the evolution of the market. However, equal weight is given to the statistics obtained at each of the even steps shown above. Whether this is a good picture of real-life trading is another matter altogether, which we shall not discuss here.

### 4.10 MINIMAL1 in action

The test of this model was to see whether it would produce better results than MINIMAL. We still used the method of averaging over `DISP_INT` time steps to obtain our statistics, as two million time step results were too many to deal with easily. However, it was possible, by a further modification to the program, to produce code that instead of averaging over `DISP_INT`, wrote the state of the system at points 2, 4, 6 and 8 in the above sequence of steps, so that each `getstats()` call was followed by the recording of the data in a file. This allowed one to produce a graph of how various quantities such as $`B`$, $`O`$ and $`Spr`$ varied with each and every time step. Besides being a satisfying thing to look at, such a graph gave a visual demonstration of the simulation at work, and an assurance that the simulation was working as it ought (Figs. 4.15 and 4.16).

The graphs show that bid and offer can meet at the same point at certain times, just before the annihilation of traders at the reaction front. Fig. 4.16 demonstrates this especially well. $`Spr`$ gradually gets closer to zero by the diffusion of the traders at the reaction front, until eventually it reaches zero and annihilation occurs, after which the best bid and best offer immediately ‘snap back’ to the position of the next best bid/offer. This is precisely the way the model is expected to work. Satisfied that our new simulation is functional, we press on and analyze results obtained thereby.

### 4.11 MINIMAL1 results

We shall perform the same analysis for MINIMAL1 results as we did for MINIMAL. To ensure the approximate validity of the approximation $`\mathrm{ln}w\mathrm{ln}cL`$, we shall take only results for which $`\mathrm{ln}w<3`$, because if we take $`c0.5`$, then for $`L=1000`$, $`\mathrm{ln}cL6`$. Values of $`\mathrm{ln}w3`$ cannot *really* be considered small compared with $`\mathrm{ln}cL`$, but if we restrict our range any further, then we would be making use of a small number of results (10 points or so) in which case reliability of results would suffer. Thus, we will use points for which $`\mathrm{ln}w<3`$, bearing in mind that the points towards the top left of the graph are less accurate.

It was decided that the last 12 points would be used to plot a graph from which the gradient and intercept would be obtained. Two factors influenced this decision—the desire to have as many points as possible to increase reliability of the results, and the need to restrict ourselves to values of $`\mathrm{ln}w<3`$—and a compromise was reached. This state of affairs was far from perfect, as the statistics generated from a mere 12 points were really not very convincing, and there were points towards the top left of the graphs for which $`\mathrm{ln}w`$ was greater than 3 or so. Nevertheless, we present the results here for comparison with later results (Figs. 4.17 and 4.18).

The graph for $`c`$ appears to suggest $`cO(0.01)`$ or thereabouts, which is much smaller than what we had before, from the MINIMAL results. The reason for this might be that 12 points do not reliably fix the gradient on a graph. Note, also, that if we take $`c0.01`$, we have $`\mathrm{ln}cL2.3`$, which requires $`\mathrm{ln}w^24.6`$. This was grossly violated in our graphs, as we allowed $`\mathrm{ln}w^2`$ to go up to 6 or slightly beyond! Thus, these results are invalid and useful information cannot be drawn from them. The statistics from these graphs actually point to their own invalidity. It would be impossible to restrict our results to $`\mathrm{ln}w^24.6`$ as our results are all $`\mathrm{ln}w^2>2`$ which already violates the regime. This forces us to take another approach.

To overcome the restriction of the range of allowed $`\mathrm{ln}w`$, we can plot $`Jw^2`$ against $`\mathrm{ln}w`$ for each value of $`L`$. This I have done, and the results are presented below:

| L | Gradient | Intercept | $`R^2`$ |
| --- | --- | --- | --- |
| 1000 | $`0.16\pm 0.02`$ | $`1.01\pm 0.13`$ | 0.559 |
| 900 | $`0.15\pm 0.04`$ | $`0.98\pm 0.13`$ | 0.348 |
| 800 | $`0.15\pm 0.04`$ | $`1.01\pm 0.12`$ | 0.383 |
| 700 | $`0.14\pm 0.10`$ | $`1.08\pm 0.28`$ | 0.0799 |
| 600 | $`0.18\pm 0.10`$ | $`1.16\pm 0.26`$ | 0.112 |
| 500 | $`0.12\pm 0.06`$ | $`0.92\pm 0.14`$ | 0.144 |

¿From Eq. (4.7), $`\mathrm{\Omega }`$ should be $`1/\pi 0.32`$. We plot the gradient as a function of $`L`$ in Fig. 4.19.

There are three sets of points: omega, min omega and max omega. The min and max series represent the minimum and maximum values omega can take, within the bounds of experimental error. Thus, they can be thought of as the ends of (imaginary) error bars. The horizontal line in the middle of the graph at $`\mathrm{\Omega }=0.32`$ represents the theoretical value. One can see that all the results agree with this predicted value, albeit with rather poor precision. Coupled with the earlier results from MINIMAL, this is mildly encouraging.

¿From the intercepts, we can work out $`c`$. We may, with some confidence, take $`\mathrm{\Omega }`$ to be $`1/\pi `$. The values of $`c`$ obtained thus were plotted against $`L`$, in Fig. 4.20.

Because of the exponentiation of the intercept involved in obtaining $`c`$, errors in the intercept were greatly magnified, leading to relatively poor precision in our determination of $`c`$. However, it agrees with the earlier result of $`c0.5`$, though this set of results seems to point, rather vaguely, to $`cO(1)`$.

### 4.12 MINIMAL2: Cornell’s approach

S. Cornell, in his paper , performed extensive simulations of diffusion-annihilation reactions, and described in detail his method of tackling the problem. His method, which allows particles of *both* species to diffuse over the same time steps, can lead to crossed-over buyers and sellers, which need then to be removed, along with the usual overlapping traders. One possible advantage of such a method, compared with MINIMAL and MINIMAL1, is that the asymmetry between buyers and sellers, in which one type of trader always diffuses before the other type, is removed, and both are allowed to diffuse simultaneously. However, an inevitable consequence of this, when used with discrete time, is the crossing over of traders. At first, one might be tempted to think that such a scenario cannot occur in the market because traders generally try to sell as high as possible, and buyers buy as low as possible. However, because exchange of information via the trading screen is not instantaneous (the screen is updated every so often), and humans can sometimes make mistakes, it is possible to have the best bid at a higher price than the best offer. In this case, usually the cross-over would be spotted by a broker, who would then inform the two traders of this situation so that they can do a deal immediately. The crossing-over characteristic of Cornell’s model, which is perfectly legitimate in chemical reactions where annihilation is not immediate but follows an exponential decay law, is actually not as far-fetched as it might seem, when applied to financial markets.

Strictly speaking, this is not the model described in , but a quick modification of the program was performed to see what such a method would yield. After several runs of this method the results did not look much different from those obtained using MINIMAL1. Further investigation is required before any conclusions may be drawn.

### 4.13 The scaling law for $`NUM`$

One interesting scaling law not mentioned in  concerns the equilibrium number of traders in the market, or $`NUM`$, where it is understood to mean the value at equilibrium. This is clearly shown in the two graphs, Figs. 4.21 and 4.22.

Fig. 4.21 shows that $`NUMJ`$ and the gradient is the coefficient of proportionality. From plotting many such graphs of $`NUM`$ against $`J`$ for different $`L`$, one can collect the gradients ($`NUM/J`$) together and thus plot $`\mathrm{ln}(NUM/J)`$ against $`\mathrm{ln}L`$ (Fig. 4.22). The graph is a straight line, with gradient $`1.999\pm 0.007`$ and intercept $`0.688\pm 0.004`$. We may infer that $`NUM/JL^2`$ with a constant of proportionality of 1/2 (from $`\mathrm{exp}(0.688)`$). Note the remarkably good correlations on the graphs and the negligible experimental errors in the determination of the gradients. We may summarize, with great confidence, that

$$NUM=\frac{JL^2}{2}$$
(4.8)
obtained empirically.

Dimensional analysis tells us that $`NUM`$ is dimensionless. The right hand side of Eq. (4.8), however, is *not* dimensionless. We can make it dimensionless by introducing $`D`$ in the denominator, as $`D/J`$ and $`L^2`$ have the same dimensions. The diffusion coefficient that we have been using is $`D=1/2`$, so the true scaling law is

$$NUM=\frac{JL^2}{4D}.$$
(4.9)
There is a firm theoretical justification for this scaling law. Consider a trader entering the system at the edge of price space. He will remain in price space until he hits the reaction front and annihilates with a trader of the opposite species. The reaction front is approximately in the middle of price space, a distance $`L/2`$ from the newly-inserted trader. From Einstein’s work on Brownian motion (see ), the mean square distance $`\lambda ^2`$ diffused by a particle in time $`t`$ is $`Dt`$. The time taken $`T`$ for the new trader to diffuse to the reaction front is given by

$$\left(\frac{L}{2}\right)^2DT.$$

Now, we multiply both sides by $`J`$, the rate of insertion of traders, and divide by $`D`$:

$$\frac{JL^2}{4D}JT.$$

$`JT`$ is simply the number of traders inserted between the time of the entrance of our trader into the market, and the time of his disappearance through annihilation. This is precisely the equilibrium number of traders in the market. One can see this in the following way: the total number of traders keeps increasing with every new trader inserted, up till a time $`T`$ after the insertion of the first trader, when annihilation starts to occur, and balances the flux of new traders. This is only approximate, of course, to within a factor of order 1. Thus, we have shown that $`NUMJL^2/4D`$.

### 4.14 Scaling law for $`\tau _S`$

Now, we turn our attention to the scaling law obeyed by the instantaneous dealing time, $`\tau _S`$, defined as the time between two consecutive annihilations. The program was modified so that it could record values of $`\tau _S`$ and calculate their average and variance. There is a vital difference between this quantity and those we have been investigating thus far: the time between sales is *not* a function of the state of the system. All the other quantities, such as bid-offer spread and midmarket, could be read off the trading screen, and they were quantities that were characteristic of the state of the market. The instantaneous dealing time, however, is different, in that at any time one does not know when the next sale is going to be, and thus what the current value of $`\tau _S`$ is. One can only record $`\tau _S`$ as a historical quantity. Furthermore, the number of readings one can take of this quantity is not fixed, but varies, depending on $`J`$ as well as the number of time steps the simulation is run for.

The program was run for the range of values of $`J`$ between 0.00004 and 0.007, which has been shown to produce $`w`$ that satisfies $`10<w<100`$, ensuring that the midmarket fluctuation, and other important lengths, are at least a factor of 10 or so away from the size of the system, or the size of the smallest unit of length. For the smallest value of $`J`$, there were only 55 values of $`\tau _S`$ obtained, whilst that number was 7140 for the largest value of $`J`$ used. Therefore, the average and variance of $`\tau _S`$ for the larger values of $`J`$ were more reliable than those for smaller $`J`$.

Originally, a graph of $`\mathrm{ln}\tau _S`$ against $`\mathrm{ln}J`$ was plotted for all the values, but while most of the points lay along a straight line, points corresponding to the 10 or so smallest values of $`J`$ were rather scattered and did not lie so well on the line. This was because of the relatively small number of data points over which the averages and variances were obtained. These points were not as reliable as the others and so it was decided to exclude them in the graph plotting, and include only those for which $`J0.0002`$, all of which had been averaged over at least 210 points. A new graph was plotted. We obtained

$$\mathrm{ln}\tau _S=(0.08\pm 0.03)(0.984\pm 0.005)\mathrm{ln}J$$
(4.10)
leading to the $`J`$ exponent being $`0.984\pm 0.005`$ and the constant in front of $`J`$ being $`1.085\pm 0.035`$. The exponent is close to the expected value of -1, but unfortunately it is not *quite* within experimental error. The constant in front is also close to the expected value of 1, but just outside experimental error. We plot a graph of the fluctuation of time to midmarket sale, and find the following:

$$\mathrm{ln}\mathrm{var}(\tau _S)=(0.34\pm 0.1)(1.97\pm 0.02)\mathrm{ln}J.$$
(4.11)
The exponent is expected to be 2, since the fluctuation, defined as the variance of $`\tau _S`$, goes as the square of time, which is inversely proportional to $`J`$, from dimensional analysis. We obtain $`0.71\pm 0.08`$ for the constant in front, for which we have no “expected” value, since intuition does not tell us how scattered the times should be. Dimensional analysis is certainly correct in telling us the exponent of $`J`$, but it cannot give us any clue about what the constant should be, apart from a possible dimensionless parameter.

It was noted that, because each of the simulations were run for a specified number of timesteps, namely 2 million, the total number of trades occurring within the time of each simulation was proportional to $`J`$. Thus, for smaller $`J`$, there were very few values of $`\tau _S`$ for small $`J`$ (as few as 55 for $`J=0.00004`$), and more for larger $`J`$ (7140 for $`J=0.007`$). The number of readings from which each average and variance is obtained should not affect the accuracy of the results, though it might lower its precision, i.e. more scatter about the “true” mean and variance. In order to overcome the problem of having a variable number of $`\tau _S`$ data points, the program was further modified so that the simulation would continue until a *fixed* number of trades had taken place. The limit was set at 1000, and the simulations terminated when 1000 trades had taken place and all their accompanying $`\tau _S`$ values were recorded. The total running time was somewhere in the region of 18 hours.

The results obtained were plotted on two graphs, one of $`\mathrm{ln}\tau _S`$ against $`\mathrm{ln}J`$ (Fig. 4.23), and the other of $`\mathrm{ln}\mathrm{var}(\tau _S)`$ against $`\mathrm{ln}J`$ (Fig. 4.24).

When plotting the graphs, it was found that the data for the smaller values of $`J`$, though lying approximately along a straight line joining them and the other points, were rather more scattered about that line than the other points were. They appeared less reliable and were therefore excluded from the two graphs presented. Only values for which $`J0.0002`$ were plotted. From the gradients and intercepts, the following information was gleaned:

$`\mathrm{ln}\tau _S`$ $`=`$ $`(0.08\pm 0.03)(0.984\pm 0.004)\mathrm{ln}J`$ (4.12)
$`\mathrm{ln}\mathrm{var}(\tau _S)`$ $`=`$ $`(0.32\pm 0.07)(1.97\pm 0.01)\mathrm{ln}J`$ (4.13)
In the case of $`\tau _S`$, the constant in front is $`1.086\pm 0.03`$, whilst the corresponding constant is $`0.73\pm 0.05`$ for $`\mathrm{var}(\tau _S)`$. One can see a remarkable similarity between these results in Eqs. (4.12) and (4.13), and the earlier ones in Eqs. (4.10) and (4.11). In fact, the only real difference is the reduction in the experimental uncertainty in the parameters. Thus, the second, more detailed, set of results not only confirms the validity of the first, but adds to it by narrowing the margins of uncertainty.

It is fair, in the light of such compelling evidence, to believe that $`\tau _S`$ goes as $`1/J`$, as both the exponent and the constant were found to be *very nearly* -1 and 1, respectively. Such a conclusion is supported also by intuition, which tells us that $`\tau _S=1/J`$. The constant in front has also been confirmed. For $`\mathrm{var}(\tau _S)`$, the exponent can be taken to be -2, while the constant in front is $`0.73\pm 0.05`$. It seems very likely that neither quantity has a logarithmic correction.

### 4.15 Scaling law for $`\tau _{\mathrm{reduced}}`$

There is another way of defining the instantaneous dealing time. It is identical to that given in §4.14 for single annihilations, but differs for multi-annihilation. We can call this $`\tau _{\mathrm{reduced}}`$, a sort of “reduced” instantaneous dealing time. If the last annihilation contained $`n`$ simultaneous trades, and the next one contained $`m`$ simultaneous trades, then

$$\tau _{\mathrm{reduced}}=\frac{\tau }{n+m1}.$$
(4.14)
$`n`$ and $`m`$ may be considered the respective degeneracies of the trades. It would be interesting to see what statistics this new quantity obeys.

We repeated the simulations for this reduced $`\tau `$, doing them for exactly the same range of $`J`$ as we did for normal $`\tau `$. 1000 data points were obtained for each $`J`$, making each point as reliable as any other. A graph of $`\mathrm{ln}\tau _{\mathrm{reduced}}`$ against $`\mathrm{ln}J`$ was plotted (Fig. 4.25).

This graph is interesting in that the points together form a precise *curve* rather than a line. There is very little scatter, except in those points with the smallest $`J`$’s. The rest of the points trace out a curve that tends to bend upwards as $`J`$ is increased. If we only take the points with the smallest $`J`$ in isolation, we obtain a gradient of –0.84. If we take the points with the highest $`J`$ in isolation, we obtain –0.58 instead. Thus, there is some variation of the exponent with $`J`$. If by brute force we calculate a line of best fit, its statistics would be

$$\mathrm{ln}\tau _{\mathrm{reduced}}=(0.793\pm 0.009)\mathrm{ln}J+(1.71\pm 0.07)$$

which gives an exponent of –0.793 and a “constant” of proportionality of $`5.5\pm 0.4`$.

Similarly, with the fluctuation (i.e. variance) in $`\tau _{\mathrm{reduced}}`$, the graph is not a straight line but a curve that curves upwards as $`J`$ is increased. This is shown in Fig. 4.26.

Its gradient varies from –0.47 at its absolute minimum to –1.73 at its maximum. The mean gradient and intercept are

$$\mathrm{ln}\mathrm{var}(\tau _{\mathrm{reduced}})=(1.32\pm 0.04)\mathrm{ln}J+(5.3\pm 0.3)$$

giving a very imprecise estimate of the constant as $`210\pm 70`$. The explanation of why the points do not fit a straight line well has not been found, and further investigations are needed before definite conclusions can be drawn. We will not pursue this question here.

### 4.16 The minimal model with asymmetric fluxes

It was now decided to generalize the situation to one in which asymmetric fluxes of buyers and sellers are permitted. Intuitively, one would expect that, in the steady-state, the midmarket moves linearly with time towards one end or the other, away from the edge with the higher flux. This is a particularly productive phenomenon to investigate, since it will tie in very well with our later work on the two-liquid model, where different fluxes of limit-order and market-order traders will be introduced, in order to reproduce the well-known widening of the bid-offer spread in the prelude to a crash. The analytical result, Eq. (3.5), for the speed of the moving midmarket, is suggested in  . This is the result we shall test.

Initially, the program was run for very long times, in order to get a feel for what happened to the midmarket when the fluxes were asymmetric. It was observed that the midmarket did move, though it did not always do so at a constant speed. To illustrate this, two graphs of midmarket variation are shown: Fig. 4.27 and 4.28.

Fig. 4.27 appears to be quite straight and linear, but Fig. 4.28 does not. The latter exhibits some irregularity. In fact, there were many other examples of such graphs, where the midmarket movement was initially quite linear, but after some time, deviated from linearity. Most tended to bend downwards, indicating a higher speed of movement. As the scaling law Eq. (3.5) applies only to small times, $`T<L^2/D`$, we should restrict our attention to times less than 2 million, and for these, the graphs were approximately linear.

It is significant that the introduction of asymmetric fluxes of buyers and sellers has not caused any change in the bid-offer spread. A graph was plotted of the variation of the best bid and best offer with time (Fig. 4.29), from which it may be seen that the spread has *not* increased with time. In addition, a graph of spread itself was plotted as a function of time (Fig. 4.30). Spread appears to stay constant with time (after steady-state was reached).

#### 4.16.1 $`\mathrm{\Delta }J`$ dependence

Our first task was to investigate the dependence of the speed of midmarket movement on $`\mathrm{\Delta }J`$. This was done by fixing all the other variables ($`J`$ and $`L`$) and allowing $`\mathrm{\Delta }J`$ alone to vary. A graph would then be plotted of $`\mathrm{ln}|dx/dt|`$ against $`\mathrm{ln}\mathrm{\Delta }J`$, from which the exponent of the $`\mathrm{\Delta }J`$ dependence could be deduced. A wide range of values for $`\mathrm{\Delta }J`$ were tried, but it was found that not all such values produced good results. For $`\mathrm{\Delta }J/J`$ that was greater than 0.4 or so, it was seen that the midmarket moved so quickly that it hit rock bottom (price = 0) in very little time, sometimes less 1 million time steps, which would have been insufficient time for the system to come to any sort of steady-state. In the past, it was observed that 1 million time steps or so were required before the system could come to equilibrium; now, although equilibrium can never be reached in this asymmetric situation, a sort of steady state can reasonably be achieved in that time. At the other end of the scale, it is possible to choose a $`\mathrm{\Delta }J`$ to $`J`$ ratio that is too small. The midmarket drifts so slowly that the gradient of the graph is negligibly small. Although not a problem in itself, when one plots the logarithm of this quantity, it is a large negative number. A small uncertainty in the gradient is magnified into a large uncertainty in the logarithm. In practice, what this means is that the points for which the ratio is smaller than 0.02 or so are quite scattered and do not provide much useful information. Therefore, we shall limit our investigations to the range $`0.02<\mathrm{\Delta }J/J<0.4`$.

Our choice of $`J`$ is still subject to the restrictions described in the earlier sections, that the lengths produced by such a $`J`$ (e.g. $`w`$, $`Spr`$) must be much larger than $`a`$, the lattice spacing, and much smaller than $`L`$. Thus, we require $`0.007<J<0.00004`$. In order to allow $`\mathrm{\Delta }J`$ to range over as large an interval as possible, it is wise to choose a large $`J`$, as $`\mathrm{\Delta }J`$ is limited by it. It was decided to use $`J=0.005`$, and $`L=1000`$, as usual. $`\mathrm{\Delta }J`$ was allowed to range from 0.0001 to 0.002. The graph is shown in Fig. 4.31.

As one can see, there is noticeably more scatter in the points towards the lower end of the graph (small $`\mathrm{\Delta }J`$). The statistics from the graph indicate that

$$\mathrm{ln}|dx/dt|=(0.98\pm 0.05)\mathrm{ln}\mathrm{\Delta }J(1.7\pm 0.4).$$

Thus, the observed exponent of $`\mathrm{\Delta }J`$ is in accordance with the theoretical value of 1. The constant of proportionality appears to be 0.18, with the minimum being 0.12 and the maximum 0.27.

There is no reason why we should not try to get more accurate and precise statistics from our results, which are perfectly valid. Looking at the graph, one can see that apart from the first 10 points or so, the data do fit a line very well, and that the fit is spoilt by the scatter in the lower 10 points. Therefore, we plot another graph, Fig. 4.32, of the same two quantities, this time without the first 10 points.

The statistics are:

$$\mathrm{ln}|dx/dt|=(0.99\pm 0.05)\mathrm{ln}\mathrm{\Delta }J(1.61\pm 0.35).$$

The exponent is now even closer to the expected value of 1. The constant in front is 0.20, with the limits being 0.14 and 0.29. In the light of such compelling evidence, we can conclude that there is no logarithmic dependence of any sort on $`\mathrm{\Delta }J`$ and that $`|dx/dt|`$ is proportional to $`\mathrm{\Delta }J`$. Moreover, we can now plot a graph of $`|dx/dt|`$ against $`\mathrm{\Delta }J`$ directly, from which we hope to obtain a more precise estimate of the constant of proportionality. This was done in Fig. 4.33.

Because we are no longer plotting the logarithm of small values of $`|dx/dt|`$, it is possible to re-include the first 10 points in our graph. The graph gives us an intercept of $`2.4\times 10^6`$ with an uncertainty of $`4.4\times 10^6`$, which makes the intercept effectively zero. The gradient is $`0.209\pm 0.005`$. Thus, plotting a direct graph can give us a *much* more precise estimate of the constant of proportionality.

#### 4.16.2 $`J`$ dependence

Our next task is to determine the $`J`$ dependence of the formula for the speed of a moving midmarket. An inverse dependence is expected. The limits on the ratio $`\mathrm{\Delta }J`$ to $`J`$ dictate that the maximum value of $`J`$ used be $`J_{max}=50\mathrm{\Delta }J`$ and the minimum be $`J_{min}=2.5\mathrm{\Delta }J`$. We require that $`J_{max}<0.007`$ in order to keep all lengths at least a factor of 10 away from $`L`$ or $`a`$, implying $`\mathrm{\Delta }J<0.00014`$. Thus, in the interests of maximizing the range of $`J`$ available to us, whilst using a “round” number, it was decided to set $`\mathrm{\Delta }J=0.0001`$, with $`J_{min}=0.00025`$ and $`J_{max}=0.005`$.

The results were plotted on a graph of $`\mathrm{ln}|dx/dt|`$ against $`\mathrm{ln}J`$ (Fig. 4.34).

This gave us the following:

$$\mathrm{ln}|dx/dt|=(1.01\pm 0.07)\mathrm{ln}J(16.1\pm 0.5)$$

from which a constant of proportionality of $`9.9\times 10^8`$ could be obtained (limits $`5.9\times 10^8`$ and $`16.5\times 10^8`$). Once again, the exponent confirms the analytic result. In order to refine our estimate of the constant of proportionality (which will be of use later), we do a direct plot of $`|dx/dt|`$ against $`1/J`$. Here, it was noticed that the 5 points corresponding to the smallest $`J`$ values we used were very scattered on that plot, so it was decided to exclude them (Fig. 4.35). It does indeed yield a more precise estimate: $`(9.3\pm 0.9)\times 10^8`$. The intercept, found to be $`(0.5\pm 1.1)\times 10^5`$, is effectively zero.

Therefore, we may conclude with confidence that the speed is inversely proportional to $`J`$.

#### 4.16.3 $`L`$ dependence

Finally, we test the $`L`$ dependence of the formula. This is, again, expected to be an inverse relationship. We tried many different values of $`L`$, from 100 up to 2000. Initially, $`J=0.001`$ and $`\mathrm{\Delta }J=0.00002`$ were used, but it was observed that the results were unusable for small $`L`$, such as 100 or 200, because there was far too much fluctuation from the mean midmarket. This is only to be expected from our earlier work on the fluctuation of the midmarket, $`w^2`$, which is known to depend on $`1/J`$. Therefore, it was decided to increase $`J`$ whilst keeping the ratio $`\mathrm{\Delta }J/J`$ the same in order to reduce the fluctuation from the mean midmarket.

$`J=0.005`$ and $`\mathrm{\Delta }J=0.0002`$ were used, which allowed the ratio $`\mathrm{\Delta }J`$ to $`J`$ to stay within the allowed limits (see §4.16.1) for $`100L2000`$. A log-log graph was plotted (Fig. 4.36).

The line is fairly straight, but with some scatter, especially for small $`L`$. The reason is that for these small $`L`$, the midmarket moved towards zero so quickly that not much time elapsed before it crashed into the lower boundary. There were not many points from which one could derive a gradient, and in any case, there was much fluctuation in the instantaneous speed of the moving midmarket, and sometimes the fluctuations dominated, leaving the drift barely observable. All these factors made it very hard to determine $`|dx/dt|`$. The gradient of the log-log plot is -1.2, which is not too bad, though we can certainly do better by excluding from our graph the first 6 points with the smallest $`L`$, which are the most scattered. This amounted to ignoring all the results for $`L<300`$. This new graph is shown in Fig. 4.37.

The improved statistics are:

$$\mathrm{ln}|dx/dt|=(1.02\pm 0.06)\mathrm{ln}L(3.0\pm 0.4).$$

The exponent is now much closer to the theoretical value of -1. The constant in front is $`0.05\pm 0.02`$. A more precise estimate may be obtained by a simple plot of $`|dx/dt|`$ against $`1/L`$, as we discovered earlier. Once again, the 6 points corresponding to $`L<300`$ (now appearing last) will be excluded from the plot, which is shown in Fig. 4.38.

The intercept is $`(0.9\pm 7)\times 10^6`$ which is basically equivalent to zero. The slope is $`0.045\pm 0.004`$, thus giving us an estimate of the constant of proportionality that is 5 times more precise than that provided by the log-log plot. We are therefore satisfied that the $`L`$ dependence really is an inverse one.

#### 4.16.4 Constant of proportionality

Having established the functional dependence of the speed of the moving midmarket, $`|dx/dt|`$, on $`J`$, $`\mathrm{\Delta }J`$, and $`L`$, we have demonstrated the truth of the formula suggested by Kogan and Eliezer, to within a constant of order 1. We shall now attempt to determine what that constant is. It is not necessary to calculate numerically the $`D`$ dependence (which we have never done for any quantity) because dimensional analysis alone tells us that if the speed, with units of dollars/sec, is proportional to $`\mathrm{\Delta }J/JL`$, there *must* also be a factor of $`D`$ there, in order to balance the units.

For each of the dependences we investigated, a partial constant of proportionality, which included in it the other two variables, was obtained. ¿From each of these we may infer an estimate of the constant in front of all three variables, i.e. the constant $`k`$, defined

$$x_0(t)=\frac{kD\mathrm{\Delta }J}{JL}t.$$

¿From the first set of results, concerning the $`\mathrm{\Delta }J`$ dependence, we deduce that $`k=2.09\pm 0.05`$, whilst for the second and third data sets, we obtain $`1.86\pm 0.18`$ and $`2.25\pm 0.2`$, respectively. We summarize this with a table:

| Set | $`k`$ | Min | Max |
| --- | --- | --- | --- |
| 1st ($`\mathrm{\Delta }J`$) | $`2.09\pm 0.05`$ | 2.04 | 2.14 |
| 2nd ($`J`$) | $`1.86\pm 0.18`$ | 1.68 | 2.04 |
| 3rd ($`L`$) | $`2.25\pm 0.20`$ | 2.05 | 2.45 |

Although not tremendously precise, the data do tell us that the constant is close to 2, which was the factor suggested in Eq. (3.37) of . Our most precise set of results was the 1st set, without a doubt, because *all* the points lay on a straight line, for both the direct and the log-log plots, and it is reasonable that we should pay more attention to its conclusions than to those of the subsequent sets, for which certain wildly inaccurate points had to be excluded. We can therefore conclude that our numerical simulations favour a constant of about 2.1, with an error of 0.1 or so.

This is supported by the theory. The formula in  was obtained from Eq. (3.36) in  by setting $`\zeta (x,t)`$, the density difference between buyers and sellers, to zero for the midmarket, and then solving the resulting quadratic in $`x`$ with the stochastic part (the Fourier series with noise-dependent coefficients) neglected. Expanding the square root in the quadratic formula gives Eq. 3.5. Of course, for the two-liquid model, when the flux of market-order traders exceeds that of limit-order traders, the system decouples into two reaction fronts, between limit and market order traders, with asymmetric fluxes. Here, the bid-offer spread widens with a speed equal to *twice* that calculated here for the minimal model, because there are two fronts. We will investigate the two-liquid model in the next section.

## 5 Computer simulations of the two-liquid model

Having spent much of our time and effort investigating the minimal model, it is now appropriate to turn our attention to the two-liquid model. It is more complicated than the minimal model to investigate, because there are four kinds of traders: limit-order buyer, limit-order seller, market-order buyer and market-order seller, and for this reason it is less analytically tractable. Many of Cardy et. al.’s methods which worked well with the minimal model cannot be applied here. Thus, there is good reason to resort to numerical simulations to solve the two-liquid model.

The intention of the two-liquid model is to provide a more realistic description of a financial market in the prelude to a crash. In particular, what motivated the inclusion of the “second liquid” was the desire to reproduce the well-known phenomenon, observed by market practitioners, of the widening of the bid-offer spread. It has been demonstrated (Fig. 4.29 and 4.30) that the asymmetric minimal model, though able to produce a steadily moving trade price, does not exhibit bid-offer spread widening. Market-order traders, believed to be the missing element, were added to reproduce this phenomenon. The rules of trading can be represented by the following equations, if we take $`B`$ and $`S`$ as the limit-order traders, and $`B^{}`$ and $`S^{}`$ as the market-order ones:

$`B+S0`$
$`B+S^{}0`$
$`B^{}+S0`$
Note that $`B^{}+S^{}0`$ does not occur, since market-order traders do not put up their prices on the trading screen and so they cannot see one another. In numerical simulation, the problem arises of which ones to annihilate when one has more than two types of traders at a reaction front. For example, if there are $`B`$ and $`B^{}`$ at the same point in price space as $`S`$, do we allow the limit-order traders to trade first, or let the trades take place randomly? In the interests of computational efficiency, it was decided to allow limit-order traders to trade first, and then let the remainding LO buyers trade with MO sellers. The three types of annihilations occur in the order shown above. Although this might not be representative of a real market, it should suffice for a first attempt at simulating one.

There are so many things in this model worthy of investigation that our research could do no more than scratch the surface, since our time was limited. Nevertheless, we did manage to get some interesting results, which are described in the following sections. We shall not dwell on the intricacies of the computer program used to perform these simulations, since the principles used here are no different from those used for the minimal model, and those we have already discussed in ample detail in §4.1.

### 5.1 Variation of spread with fraction of market-order traders

It was decided to measure the bid-offer spread as a function of the fraction of market-order traders. Because we are using a model in which only the flux of traders is controllable, and the total number of traders of either type fluctuates, we can only measure the fraction of market-order traders with respect to the relative fluxes of the two types of traders. In other words, the fraction of market-order traders (MO) is here defined as the MO flux divided by the *total* flux of both types of traders. Using this definition, it is possible to change the fraction of MO traders by changing the MO flux, while keeping the total flux the same. It will be seen that this definition of MO fraction is supported by theory (which does *not* support a definition based on trader number).

As time progresses, we know that in cases where MO flux exceeds LO (limit-order) flux, the bid-offer spread widens. The problem of what value to take as the spread for a particular flux configuration thereby arises. It was decided to define the spread here as the average spread between $`t=300,000`$ and $`t=330,000`$. This interval was chosen because it was *believed* to be sufficient time for the system to come to a sort of steady-state, yet not too long so that the two separating reaction fronts crash into the boundaries of price space. The mean flux was chosen to be 0.0005, a relatively high flux, because midmarket fluctuations, $`w^2`$, which we wanted to minimize in order to have more accurate results, are known to go as $`1/J`$. Thus, a 10% MO flux would correspond to $`J_{MO}=0.0001`$ and $`J_{LO}=0.0009`$. Simulations were performed for many different fractions and the results were plotted on a graph (Fig. 5.1).

The critical point seems to be approximately 0.5, above which the spread begins to grow very rapidly Preliminary results that bifurcation point in the Two-Liquid Model occurs at $`f=1/2`$ was obtained by Adrian McGowan and one of us (I.I.K.) in the spring 1999.. There is some fluctuation in the results, but they do fit a curve quite well, though they would fit a straight line probably just as well—we cannot tell at this stage. The value of the spread is undefined for 100% market-order traders.

It is interesting to consider what would happen to the spread as the MO fraction tends to 1. We may apply Eq. 3.5 to the two-liquid model if we assume that it consists of two decoupled minimal systems, of MO against LO traders, with reaction fronts which move apart (this only works if $`J_{MO}>J_{LO}`$). We must also assume that the two fronts act independently of one another. This is only an approximation and does not apply at the beginning of the simulation, when the system is still trying to reach a steady-state. Eq. 3.5 tells us that for a system with asymmetric fluxes $`\overline{J}+\mathrm{\Delta }J`$ and $`\overline{J}\mathrm{\Delta }J`$, the speed of movement of the midmarket (which is effectively the same as the rate of movement of the reaction *front*, since the spread stays constant) is $`2D\mathrm{\Delta }J/\overline{J}L`$. In our case, $`J_{MO}=\overline{J}+\mathrm{\Delta }J`$ and $`J_{LO}=\overline{J}\mathrm{\Delta }J`$. Making the necessary substitutions, and remembering that the fraction $`f`$ of MO traders is defined $`f=J_{MO}/(J_{MO}+J_{LO})`$, we obtain

$$\frac{d(Spr)}{dt}=\frac{4D(2f1)}{L}$$

which tells us how the spread increases with time. Integrating, we have

$$Spr(t)=\frac{4D(2f1)t}{L}+S_0$$
(5.1)
where $`S_0`$ is the “natural” spread which exists in the absence of the drift, at $`f=0.5`$. The formula predicts, therefore, that the graph of spread against fraction $`f`$ would be linear for $`f>0.5`$, with a gradient of $`8Dt/L`$. In our case, this would have been 1260, if we use $`t=315,000`$. The spread at $`f=0.5`$, $`S_0`$, is about 84, leading to a prediction of 714 for the y-intercept.

We have already observed that the portion of the graph (Fig. 5.1) for $`f>0.5`$ is more like a curve than a straight line. However, if we use a straight line as a first approximation, and fit it to the points $`f>0.5`$, we can see how the theoretical prediction matches up with the data. This was done in Fig. 5.2 which was a plot of all the points for which $`f0.5`$, excluding the last three points which were obvious outliers.

The graph gives a gradient of $`1354\pm 64`$ and an intercept of $`654\pm 48`$. Comparing the gradient with the theoretical value of 1260, one can see that the agreement is quite good. The y-intercept, too, is in surprisingly good agreement with the calculated value. Having said that, it is clear that Eq. (5.1) gives little more than an estimate of what the spread is. In any case, the equation only works in the region $`f0.5`$.

A market practitioner would want a better fit to the data than a crude straight line. To this end, we have fitted the curve (less the 3 outliers) to a cubic, as shown in Fig. 5.3. The fit is remarkably good.

It cannot be determined with any certainty at this stage whether the portion of $`f>0.5`$ is a curve, like a cubic, or a straight line.

One wonders whether, in spite of the excellent cubic fit performed above, a single description of the variation of spread, or indeed any quantity, can be valid for all $`0<f<1`$. A cause for concern lies in the very sharp transition from a stable system $`f<0.5`$ to an unstable one $`f>0.5`$, with $`f=0.5`$ being the critical point. We expect 0.5 to be the critical point, because an infinitesimal increase of $`f`$ beyond this point would produce an imbalance of MO and LO fluxes, which would, in our model, inevitably lead to a steady widening of the spread. Yet is this the case? We see a sharp rise in spread at around 0.5, but does a phase transition really occur at exactly 0.5? The spread plotted was only the average around $`t=315,000`$; perhaps we have been observing the spread too soon after the beginning of the simulation, and the system has not had time to reach a steady-state? These are all valid questions, which we now try to answer by performing more extensive simulations.

### 5.2 More definite results for spread

Admittedly, $`t=300,000`$ to $`t=315,000`$ is rather soon to start observing market behaviour; we know that for symmetric fluxes, the market requires approximately 1 million time steps to come into equilibrium. We were forced to do so in order to maximize the range of $`f`$’s accessible to us. It was necessary to compare spreads measured at the *same point in time* for each simulation. Had we chosen a later time to observe, e.g. 1 million steps, we would not have been able to observe the variation of spread with $`f`$ for $`f0.75`$ or so, for which the spread at $`t=10^6`$ would have exceeded 1000One might think that the problem can be solved by simply increasing $`L`$. Not only does this increase the size of price space, but it also slows down the widening process (speed of widening $`1/L`$) and overall, the time taken for the spread to reach $`L`$ goes as $`L^2`$. However, the time taken for the system to reach equilibrium as a whole goes as $`L^2`$, cancelling exactly with the increase in time gained from increasing $`L`$.. If we look at $`f=0.0`$, and compare the spread obtained by averaging from $`t=300,000`$ to $`t=330,000`$, and that obtained by averaging from $`t=10^6`$ to $`t=2\times 10^6`$, we see that they are 45.5 and 41.7, respectively. An acceptable difference, if we bear in mind the inherent fluctuations. However, for $`f=0.50`$, the spreads are 84.0 and 127.5, respectively. They differ by a factor of 1.5. Therefore, it is probable that the system had *not* reached a steady-state by the 300,000th time step.

In order to probe this further, the spread calculations performed earlier were repeated, this time taking the spread as the average spread between 1 million and 2 million time steps (Fig. 5.4). This is the most comprehensive set of results for two-liquid spread so far.

Looking at the graph, one is immediately struck by how well the points $`f<0.5`$ fit a straight line. The reason for this is that the average over 1 million time steps is less prone to fluctuation than that over a mere 30,000 steps. We are limited in what region of $`f`$ we can probe, since the spread shoots upwards very quickly once $`f`$ exceeds 0.5, and reaches 1000 around 0.64 or so. One can also see that the graph divides neatly into two straight sections: $`f<0.5`$ and $`f>0.5`$. This is encouraging, as it accords with what we predicted using the formula for asymmetric fluxes applied to the two-liquid model, Eq. (5.1). Once again, we plot the two sections, $`f<0.5`$ and $`f>0.5`$ separately to determine the gradient and intercept (Fig. 5.5 and 5.6).

We consider first the section corresponding to $`f0.3`$, which is the straightest part of the $`f<0.5`$ half. The full plot in Fig. 5.4 shows an almost perfect straight line. The expanded version in Fig. 5.5 supports this, with the fitted line achieving a correlation coefficient in excess of 0.94. The gradient and intercept are, respectively, $`47.6\pm 2.2`$ and $`40.4\pm 0.4`$. Initially, attempts were made to include points up to $`f=0.4`$ but it was found that the points did not lie so well on a straight line. At the moment, we do not have a theory that can explain this straight portion. Note that the spread at 0.5 cannot be worked out by extrapolation from this formula, since the spread starts to deviate from this line at $`f0.4`$.

Turning to the $`f0.525`$ portion, we note first that, contrary to previous results, the points quite clearly fit a line. The gradient and intercept are $`5927\pm 223`$ and $`2928\pm 130`$, respectively. Theory predicts that, for $`t1,500,000`$, the gradient $`d(Spr)/df`$ would be $`8Dt/L=6000`$. The intercept is $`S_04Dt/L`$, the first term of which has to experimentally determined. We find that the equation describing $`f0.5`$ is

$$Spr(f)=(5927\pm 223)f(2928\pm 130).$$

The gradient agrees with theory, to well within experimental error. From Fig. 5.4, we find, by visual inspection, that the spread at the point of symmetry, $`f=0.5`$, is $`140\pm 20`$. Incidentally, *neither* straight line section, when extrapolated to $`f=0.5`$, gives the correct result, since 0.5 is in the curved section joining the two straight sections, and is larger than either straight line formula would predict. Now that we know $`S_0`$, we also know what theory predicts for the intercept of the $`f0.525`$ graph: it is $`2860\pm 20`$ (the uncertainty comes from our experimental determination of $`S_0`$). The actual intercept agrees well with our prediction. Thus, we are pleased to see good agreement with theory.

In the light of these much more detailed results, especially Fig. 5.4, it appears that the spread starts to rise rapidly *just before* $`f=0.5`$, because we have already observed that the spread at that point of symmetry deviates from either linear regime. We recall that, in our method of simulation, whenever there was an abundance of traders that needed to be annihilated at a point in price space, we decided to always allow the LO traders to annihilate first. Such a technique would, in general, cause the best LO bid/offer to be further apart than it would be if we had not made that imposition and instead had allowed the LO and MO traders to trade with equal probability. This is because the spread is determined by LO traders alone, and if LO traders are more likely than MO traders to get annihilated, it is reasonable that the best bids and offers would ‘snap back’ more often, resulting in an increased average spread. Immediately, we can see that this preference for LO-LO trades only affects those trades where both LO buyers and LO sellers are present, i.e. when the spread is zero. Therefore, for $`f>0.5`$, where the spread widens at a constant rate, and apart from at the beginning, the LO traders are completely separated, and we would expect the LO-LO preference to make no difference, since the LO traders on either side are no longer in contact. For $`f`$ just below 0.5, we expect no steady drift of the best bid or offer due to asymmetric fluxes, and so the situation is less clear. During all those times when the best LO bid and offer are not in contact, the priority given to LO-LO trades makes no difference. When they are coincident, however, *more* of the LO traders are going to be annihilated than they would otherwise. As a result, in such an annihilation where traders of all types are present, the number of LO traders will be more greatly diminished than that of MO traders. To consider what happens to the LO best bid, envisage a situation where, at the midmarket, the number of LO buyers is less than the total number of LO/MO sellers, which in turn is less than the total number of LO/MO buyers. Following our scheme, LO-LO trades occur first, followed by LO-MO trades. All LO buyers are eliminated. The only way for the best bid *not* to change would be for the number of LO buyers to exceed the total number of sellers. Had we used a random scheme of mutual annihilation, however, since buyers outnumber sellers, there would be no reason why some LO buyers should not be left. Similarly, if we perform this analysis for sellers, we would find a similar prejudice against the best offer staying where it is. In both cases, the number of LO traders has to outnumber the total number of the opposite type of trader in order for the best bid/offer to remain unchanged. For a random annihilation scheme, the only condition is that one type of trader outnumber the other, and then there is the *possibility* of either the best bid or best offer to remain unaltered. Therefore, it can be seen that increased spread is an artefact of the our simulation.

By itself, the above analysis does not explain why the critical point might occur below 0.5. To do that would require one to know quantitatively how the LO-LO preference affects the spread. Alternatively, to obviate this problem, it might be useful to change the simulation algorithm so that trades are indeed random, and that no preference is given to any type of trade. Such a modification would greatly increase processing time though it would probably resolve our uncertainty.

It is interesting to note the numerical prefactor that is defined as the ratio of the spread at the point of symmetry ($`f=0.5`$) to that for the minimal model ($`f=0.0`$). From the results we already have, we can immediately give a value for it. The former is $`140\pm 20`$ and the latter $`40.4\pm 0.4`$, both coming from the graph of spread in the region $`f<0.5`$. The ratio is therefore $`3.47\pm 0.50`$.

## 6 Conclusion

The numerical analysis performed in this paper confirmed that $`w^2`$, $`Spr^2`$ and $`\mathrm{var}(Spr)`$ are all proportional to $`D/J`$, as predicted by dimensional analysis. The results also confirmed the existence of logarithmic correction for $`w^2`$, but not for $`Spr`$ or $`\mathrm{var}(Spr)`$. This is not surprising, because the midmarket is ergodic and, given enough time, it will cover the entire width of price space. Increasing $`L`$ increases the range in price space over which the midmarket can roam; thus, the midmarket variance diverges as $`L\mathrm{}`$. The constant in front of the $`D/J`$ in the $`w^2`$ scaling law was confirmed to be $`1/\pi `$. The logarithmic parameter for $`w^2`$ is $`c0.56`$, $`0.32<c<0.97`$, for $`L=1000`$. It appears that $`cO(1)`$ for $`500<L<1000`$. $`Spr`$ and $`\mathrm{var}(Spr)`$ do not seem to have logarithmic corrections. This is because $`Spr`$ (and $`\mathrm{var}(Spr)`$, which is derived from $`Spr`$) is a fundamentally different quantity from the midmarket variance. Intuitively, changing $`L`$ should not change $`Spr`$, since it is a property of the reaction front. The constants in front of $`D/J`$ for $`Spr`$ and $`\mathrm{var}(Spr)`$ are $`3.615\pm 0.065`$ and $`0.95\pm 0.02`$, respectively, whilst the approximate ‘constant’ in front of $`D/J`$ for $`w^2`$ is $`1.9\pm 0.2`$. We can therefore observe a hierarchy in the three major lengths: $`\mathrm{var}(Spr)<w^2<Spr`$ according to $`0.95<1.9<3.6`$. It was discovered that $`NUM`$ scaled as $`JL^2/4D`$, and a theoretical justification was given. This law was the most exact of all the scaling laws that the data supported; there was almost perfect correlation for all graphs of $`NUM`$ against $`J`$ and $`L`$. Time to midmarket sale ($`\tau _S`$) was investigated, and found to be equal to $`1/J`$, as expected. Its fluctuation also scaled as $`1/J`$, but with $`0.73\pm 0.05`$ as the constant of proportionality. Neither quantity showed signs of having logarithmic corrections. A further quantity, $`\tau _{reduced}`$, a sort of scaled time, was investigated. Its scaling law seemed a little less straightforward, as the points on the log-log plot quite clearly followed a curve rather than a line. Further investigation is needed to determine its exact nature. The minimal model was extended by allowing asymmetric fluxes of buyers and sellers. Simulations verified the $`\mathrm{\Delta }J`$, $`J`$ and $`L`$ dependences, and the constant of proportionality, of the formula (Eq. (3.5) in this paper) predicted by Kogan et. al. in  for the speed of a moving midmarket in the presence of asymmetric fluxes.

Investigations into the two-liquid model produced widening of the bid-offer spread in the prelude to a crash, when the flux of market order traders exceeded the flux of limit order traders. Application of Eq. (3.5) to the spread yielded approximate, though uncertain, agreement. A cubic fit to the spread as a function of MO fraction was performed. Further, more detailed, results confirmed that the graph of $`Spr`$ against $`f`$ was made up of two straight sections, in agreement with theory. The predicted gradient and intercept of the second straight section were corroborated by experiment. The formula for asymmetric fluxes was found to describe the bifurcated regime well, but sheds no light on the equilibrium regime, where $`f<0.5`$. The critical point was found to be approximately $`f=0.5`$, or perhaps just before. The ratio of spread for $`f=0.0`$ and $`f=0.5`$ was calculated to be $`3.47\pm 0.50`$.

The analysis presented in this paper has covered five major scaling laws described in . However, there are many more to be tested, and so there is much scope for further research. Other scaling laws that might be investigated include density near the best bid/offer as a function of deal rate, and higher correlation functions describing equilibrium. We have only scratched the surface of the two-liquid model; it would be interesting to see how the other scaling laws are modified in the prelude to a crash. Finally,  describes a third model, the bias model, which attempts to incorporate herd, or crowd, effects. It does so by modifying the diffusion operators to contain a drift, which depends on whether the last trade has moved upwards or downwards. This modification models momentum trading, where the traders have a tendency to go with the crowd. Such a model is designed to develop instabilities, and the drift parameter may be considered to be in a meta-stable basin, buffeted back and forth by random market movements (diffusion), until the drift is so large that it is knocked out of its basin, into the surrounding unstable region. Once there, it gathers momentum by positive feedback, leading to a crash. In future work we shall investigate the bias model.

## 7 Acknowledgements

The work of one of us (D.L.C.C.) was partly supported by the summer project grant from Physics Department, Oxford University. D.L.C.C. would like to thank Ming Lee and Stella Chan for their support and understanding. D.E. would like to thank Matt Moen, James Lee, Monique Elwell. I.I.K. would like to thank John Cardy, Juan Garrahan, Zoltan Rácz, David Sherrington and Oleg Zaboronski and acknowledge stimulating atmosphere of ESF programme SPHINX.

## Appendix A Appendix

### A.1 The master equation and the diffusion coefficient

From the master equation for the trader diffusion process, we may derive the diffusion equation and thus an expression for the diffusion coefficient $`D`$.

Let $`P(x|t)`$ be the probability of having a particle at point $`x`$. Let $`p`$ be the probability for a particle to hop *away* from a given point, in each time quantum $`\tau `$. Each hop is to a point a distance $`a`$ away from the original point of the particle. Let the probability of hopping to the left and to the right be equal, i.e. $`p/2`$. Therefore, the probability of moving to the left is $`p/2`$, moving to the right $`p/2`$, and staying put $`(1p)`$. We write the master equation by considering the probability for the particle to be at $`x`$ at a time $`t+\tau `$:

$`P(x|t+\tau )`$ $`=`$ $`{\displaystyle \frac{p}{2}}P(x+a|t)+{\displaystyle \frac{p}{2}}P(xa|t)+(1p)P(x|t)`$
$`P(x|t+\tau )P(x|t)`$ $`=`$ $`{\displaystyle \frac{p}{2}}[P(x+a|t)+P(xa|t)2P(x|t)]`$
We expand either side in a Taylor series about $`t`$ and $`x`$. The odd derivatives cancel out on the right hand side, leaving just the even derivatives, which add. For small $`a`$ and $`\tau `$, we may neglect terms of second order or higher in $`\tau `$ on the LHS and fourth order or higher in $`a`$ on the RHS. The partial derivatives are all evaluated at $`x`$ and $`t`$.

$`\tau {\displaystyle \frac{P(x|t)}{t}}+\mathrm{}`$ $`=`$ $`{\displaystyle \frac{p}{2}}\left[{\displaystyle \frac{^2P(x|t)}{x^2}}\left({\displaystyle \frac{a^2}{2}}\right)+\mathrm{}+{\displaystyle \frac{^2P(x|t)}{x^2}}\left({\displaystyle \frac{a^2}{2}}\right)+\mathrm{}\right]`$
$`{\displaystyle \frac{P(x|t)}{t}}`$ $`=`$ $`\left({\displaystyle \frac{a^2p}{2\tau }}\right){\displaystyle \frac{^2P(x|t)}{x^2}}.`$
This is the diffusion equation for $`P(x|t)`$, with diffusion coefficient $`D=a^2p/(2\tau )`$. In our simulations, $`p=1`$ so $`D=a^2/(2\tau )`$, or $`D=1/2`$, for $`a=\tau =1`$.

### A.2 Proof of $`X=\overline{X}`$ and $`\mathrm{var}(X)=\overline{X^2}\overline{X}^2`$

If we have a set of $`N`$ results $`\{X_i\}`$ and we choose to record only the averages over $`n`$ steps, $`\overline{X}_j`$, of the $`j^{\mathrm{th}}`$ set of $`n`$ values, such that

$`\overline{X}_j={\displaystyle \frac{1}{n}}{\displaystyle \underset{i=n(j1)+1}{\overset{nj}{}}}X_i,`$ $`1jN/n,`$
then we may obtain the average ($`\mathrm{}`$) over the entire set of data by

$$X=\frac{1}{N}\underset{i}{}X_i=\frac{1}{N}\underset{j=1}{\overset{N/n}{}}n\overline{X}_j.$$

This we may rewrite as

$$X=\frac{1}{N/n}\underset{j=1}{\overset{N/n}{}}\overline{X}_j=\overline{X}$$

which is simply the mean of the set of recorded averages ($`n`$ is a constant). Similarly, we may obtain the variance of an entire set of data from knowing the sum of the squares, and the mean, of sets of $`n`$ values:

$`\mathrm{var}(X)`$ $`=`$ $`{\displaystyle \frac{_iX_i^2}{N}}\left({\displaystyle \frac{_iX_i}{N}}\right)^2`$
$`=`$ $`{\displaystyle \frac{_jn\overline{X^2}_j}{N}}\left({\displaystyle \frac{_jn\overline{X}_j}{N}}\right)^2`$
$`=`$ $`{\displaystyle \frac{_j\overline{X^2}_j}{N/n}}\left({\displaystyle \frac{_j\overline{X}_j}{N/n}}\right)^2`$
$`=`$ $`\overline{X^2}\overline{X}^2.`$
Thus, the variance of an entire set of data is equal to the mean of the averaged squares minus the square of the mean of the recorded averages. Armed with these relations, we may reconstruct the mean and variance of the whole set of data from the averaged recorded data.