1 Introduction↩︎

Financial markets operate through intricately coupled collections of order books. Order books can be differentiated at the most coarse level into either lit or latent order books, both play pivotal roles in the dynamics of price emergence. Lit order books provide an environment where all market participants can readily view and engage with displayed orders [1] and actively transact. This transparency not only fosters liquidity but also enables traders to promptly identify available orders, thereby contributing to the creation of a robust and efficient market.

A distinctive attribute of lit order books is their immediate and transparent execution mechanism. Trades within lit order books transpire swiftly at the best available price, aligning with the price-time priority principle [1]. In contrast, latent order books [2] encompass hidden or non-displayed orders, providing a realm of confidentiality for traders, but really represent potential orders and demand provided by the largest long-term investors in a market.

[3]introduced the idea of using a reaction-diffusion model of the latent order book. [4] implemented a numerical solution using a stochastic finite difference method [5], [6] to simulate and then calibrate a reaction-diffusion market model. The model was then extended by [7] to consider anomalous diffusion with non-uniform sampling times. This model was then used to numerical demonstrate that the Epps effect can be generated as an emergent correlation [8]. We use this model formulation to demonstrate how correlation can emerge from two coupled lit order books and thus verify if the model can naturally generate an Epps effect as an emergent property [9]. The Epps effect has been observed in various empirical studies of financial markets [10].

Here correlation is not imposed exogenously through either correlated noise or latent multivariate diffusions. It emerges endogenously from an order–flow coupling between books. Here the coupling represents a pair–trader mechanism that responds to transient mispricing by placing side–selective orders near the reaction boundary. In contrast to Hawkes–process or multivariate diffusion models, where correlation parameters are fixed primitives [11], [12]; here the effective correlation is an output of the dynamics. The long–run correlation appears only after sufficient aggregation, while short–horizon correlation is suppressed because the coupling response is neither instantaneous nor perfectly synchronous.

In summary: we extend the numerical simulation and framework of [8] to an analytic derivation and approximation. This aims to confirm analytically that when there are interactions mediated by pairs traders who react to temporary deviations between prices and submits orders that induce mean reversion correlation naturally emerge and reproduce an Epps affect. Our goal is threefold: first, to clarify the discrete-to-continuum construction underlying the coupled equations; second, to show analytically how this interaction generates scale-dependent correlations and an Epps effect; third, to demonstrate the relationship between time-subordination from events to model time and how this relates to the affect generated by the trader couplings.

2 The coupled limit order books↩︎

We start with a discrete price grid with spacing \(\Delta x\) and a sequence of event times \(0<t_1<t_2<\dots\). At each event time, orders may be created or cancelled, and existing orders may shift by one price level. For a single book, the microscopic dynamics can be represented by a discrete-time random walk in price space with killing and source terms.

Let \(\varphi_i^{n(j)}\) denote the signed order density at grid point \(x_i=i\Delta x\) for asset \(j\) at event index \(n\). Transition probabilities \(\tfrac{1}{2}r\) for left and right moves, cancellation probabilities, and source terms define a master equation for \(\varphi_i^{n(j)}\) [7], [8]:

Here \(K_{n-m}\) is a memory kernel for a process with Sibuya waiting times [6] which defines the memory due to the fractional diffusion. There is a random driving force \(F^{(j)}_{n-1}\) and an incremental source term \(c^{(j,k)}_{i,n-1}\). The coupling between the book is a cross-book reaction term.

2.1 Coupled update equations↩︎

If there is an pair-wise order book coupling from the \(j^{th}\) order book to (say) the \(k^{th}\) order book this will be carried by the coupling term \(\ell^{(j,k)}\). If there is no coupling then we can drop the \(k\) index on the source term. Boundary condition considerations (to ensure that that the integral under \(\varphi^{(i)}\) is constant) make it convenient that we will only consider the lit order book sources [4], [7], [8]: \[\label{eq:source-term} s^{(j)}(x,t) = - \lambda_j \mu_j (x - p^{(j)}(t)) e^{\mu_j\left(x - p^{(j)}(t)\right)^2}. \\\tag{2}\] Here, using pair-wise couplings between the \(j^{th}\) and \(k^{th}\) order books [13]: \[\begin{align} \ell^{(j,k)}(x,t)=G(x,t,p^{(j)},\Delta p_{jk}) = G_j(x,t) \\ \ell^{(k,j)}(x,t)=G(x,t,p^{(k)},\Delta p_{kj}) = G_k(x,t). \label{eq:ell-crosscoupling} \end{align}\tag{3}\] The coupling equations are a function of the difference in mid-prices of two order books where, if \(p^{(j)}\) (in \(\varphi^{(j)}\)) is above \(p^{(k)}\) (in \(\varphi^{(k)}\)), then more bids are placed above the mid-price in \(\varphi^{(j)}\) to push the mid-price down to \(\varphi^{(k)}\)’s mid-price. To achieve this we define: \[g_j(x,t) = - \lambda_j \mu_j x e^{\mu_j x^2}.\] Here \(\lambda_j\) and \(\mu_j\) are constants specific to the \(j^{th}\) order book, and \(\Delta p_{jk}\) is the difference between the mid-prices of the two order books: \(\Delta p_{jk} = p^{(j)} - p^{(k)}\): \[\begin{align} G_j = \left\{ \begin{array}{ll} g_j\left(x-p^{(j)}(t)\right)\, \Delta p_{jk} &,x > p^{(j)}(t), \Delta p_{jk} > 0\\ g_j\left(\frac{j}{\Delta p_{jk}}(x-p^{(j)}(t))\right) &,x \leq p^{(j)}(t), \Delta p_{jk} > 0 \\ g_j\left(x-p^{(j)}(t)\right)\, \Delta p_{jk} &,x \leq p^{(j)}(t), \Delta p_{jk} \leq 0\\ g_j\left(\frac{1}{\Delta p_{jk}}(x-p^{(j)}(t))\right) &,x > p^{(j)}(t), \Delta p_{jk} \leq 0. \end{array} \right. \label{eq:G-coupling} \end{align} \nonumber\tag{4}\] This has the interpretation of an external agent (such as a pairs trader) observing the system, and buying (or selling) one asset according to whether the mid-price of the other asset is above (or below) some price threshold. In this way, the system can be generalised to many assets being traded by pair traders. Here we will focus on only two assets.

2.2 Moving reaction boundary↩︎

The trade prices for the \(j^{th}\) order book are: \[p^{(j)}(t) = {{x} \,\, \operatorname{s.t.}\,\, \varphi^{(j)}(x,t)=0}.\] We assume that the annihilation rates are constants: \[a^{(j)}(x,t)=\nu_j.\] The creation term is separated further into a source terms \(s^{(j)}\), coupling terms \(\ell^{(j,k)}\), and shocks \(\delta^{(j)}\) for the \(j^{th}\) assets order book: \[c^{(j,k)}(x,t)=s^{(j)}(x,t)+\ell^{(j,k)}(x,t)+ \delta^{(j)}(x,t).\] Order book shocks \(\delta^{(j)}\) are used to estimate the price impact [7]. These involve introducing a shock of size \(Q\) and measuring the change in price as a result. The source terms \(s^{(j)}\) are assumed to be either latent order book sources or lit order book sources. Lit order book source terms have vanishing boundary conditions, while the latent order book have finite boundary conditions values.

3 From discrete order flow to a continuum limit↩︎

We assume that the diffusion limit exists \[\label{eq:d} D_{\alpha} = \lim_{\substack{\Delta x \rightarrow 0 \\ \Delta t_m \rightarrow 0}}\frac{r}{2}\frac{\Delta x^2}{\Delta t_m^{\alpha}},\tag{5}\] where \(D_{\alpha}\) is some diffusion parameter. We can then use this to set the lattice price grid increments in terms of the lattice time increments following [7]. We will treat this parameter as the same across all the order books, similarly for the fractional time parameter \(\alpha\) with \(0<\alpha\le 1\), the master equation converges to a reaction–diffusion equation. For \(\alpha=1\) one obtains ordinary diffusion; for \(\alpha<1\) the limit is anomalous diffusion with a fractional time derivative, corresponding to heavy-tailed waiting times between events.

Then with some approximation function \(\varphi^{(j)}_{\Delta}(x,t)\) for the order book densities on some background lattice [7] so that \(\varphi^{i(j)}_n = \varphi^{(j)}_{\Delta}(x_i,t_n)\) with non-uniformly sampled time increments \(\Delta t_n\). The order density \(\varphi^{(j)}(x,t)\) in the continuum operational time diffusion limit satisfies \[\partial_t \varphi^{(j)} = D_\alpha \, \partial_x^2 \bigl[ D_t^{1-\alpha} \varphi^{(j)} \bigr] - \nu_j \varphi^{(j)} + c^{(j,k)}(x,t),\] where \(c^{(j)}\) collects order creation and interaction terms.

The transaction price is not imposed a priori. Instead it is defined endogenously as the reaction boundary \[p_j(t) := \{ x : \varphi^{(j)}(x,t)=0 \}.\] This moving boundary separates buy and sell dominated regions of the order book. In the continuum description, the price dynamics is therefore encoded implicitly through the evolution of \(\varphi^{(j)}\). A key step in the analysis is a local linearisation of this zero-crossing condition.

3.1 Regularised analytic coupling notation↩︎

The piecewise coupling in Equation 4 is the coupling used in the numerical simulations of [8]. For the analytic derivations below we use a regularised local-response representation of the same pair-trader mechanism. This preserves the directed, side-selective mean-reverting order flow, but removes the inverse-spread singularity at the symmetric state around which the analytic approximation is constructed.¹ Here we also change the notation so that the differences are explicit.

Let the price given by the reaction boundary represent the transaction price at continuum operational time \(t\) for the \(j\)-th order book: \[p_j(t):=p^{(j)}(t),\qquad y_j(x,t):=x-p_j(t).\] Now define the directed spread \[z_{jk}(t):=p_j(t)-p_k(t), \qquad z_{kj}(t)=-z_{jk}(t).\] The analytic source-shape kernel is written \[q_j(y)=-\lambda_j\mu_j ye^{-\mu_j y^2}, \label{eq:qj-analytic}\tag{6}\] where \(\lambda_j>0\), and \(\mu_j>0\). The derivation only requires \(q_j\) to be odd and locally linear near the mid-price. ²

The regularised pair-trader coupling from book \(k\) into book \(j\) is \[\ell^{(j,k)}(x,t) =\gamma_{jk}z_{jk}(t)q_j\!\left(y_j(x,t)\right) W\!\left(y_j(x,t),z_{jk}(t);\varepsilon\right), \label{eq:regularised-coupling}\tag{7}\] where \(\gamma_{jk}\geq0\) is a coupling-strength parameter and \(W\) is a bounded side-selection function. A smooth choice [14], [15] is \[W(y,z;\varepsilon) =\frac{1}{2}\left[1+\tanh\!\left(\frac{yz}{\varepsilon}\right)\right], \qquad \varepsilon>0. \label{eq:side-selection-function}\tag{8}\] Thus \(z_{jk}>0\) selects the \(y_j>0\) side of book \(j\), while \(z_{jk}<0\) selects the \(y_j<0\) side. The corresponding hard-threshold limit is \(W_0(y,z)=\mathbf{1}_{\{yz>0\}}\). This notation makes the spread direction, local price coordinate and coupling strength explicit, and provides the differentiability needed for the price-level approximation.

4 Analytic decomposition of the Epps effect↩︎

Let \(R^{(j)}_\Delta(t)=p_j(t+\Delta)-p_j(t)\) denote returns over scale \(\Delta\). The realised correlation [16] is \[\rho_\Delta=\frac{\mathbb{C}\mathrm{ov}(R^{(1)}_\Delta,R^{(2)}_\Delta)}{ \sqrt{Var(R^{(1)}_\Delta)Var(R^{(2)}_\Delta)}}. \label{eq:realised-corr}\tag{9}\]

Equation 9 is the population analogue of the realised-correlation estimator used in the high-frequency realised-covariation literature [16]. When the two assets are observed on non-synchronous event clocks, the relevant overlap structure is closely related to non-synchronous covariance estimation [17].

4.1 Subordination-only mechanism↩︎

If prices evolve in operational time but are observed through independent clocks, the correlation is reduced by limited clock overlap. One finds \[\rho_\Delta^{\mathrm{sub}}=\rho_\infty \, \mathcal{A}_{12}(\Delta),\] where \(\mathcal{A}_{12}(\Delta)\) is an overlap factor, and \(\rho_{\infty}\) is the instantaneous correlation of the two Brownian motions in operational time. For Poisson refresh rates this yields \[\rho_\Delta^{\mathrm{sub}} \simeq \rho_\infty \left[1-\frac{1-e^{-\lambda_{12}\Delta}}{\lambda_{12}\Delta}\right].\]

See 8 where for independent Poisson refresh clocks with rates \(\lambda_1\) and \(\lambda_2\), the pooled refresh process is Poisson with rate \(\lambda_{12}=\lambda_1+\lambda_2\). Here the waiting time until either book refreshes is the minimum of two independent exponential waiting times and is therefore exponential with rate \(\lambda_1+\lambda_2\).

4.2 Coupling-only mechanism↩︎

Linearising the zero-crossing condition gives a price-level approximation: \[\begin{align} dp_1=\sigma_1 dW_1-\kappa_1(p_1-p_2)dt, \\ dp_2=\sigma_2 dW_2+\kappa_2(p_1-p_2)dt. \label{eq:coupling-only} \end{align}\tag{10}\] The spread relaxes with rate \(\kappa=\kappa_1+\kappa_2\). Integrating the resulting cross-covariance kernel gives \[\rho_\Delta^{\mathrm{coup}} \simeq \rho_\infty \left[1-\frac{1-e^{-\kappa\Delta}}{\kappa\Delta}\right].\]

See 9 where \(\rho_\infty\) denotes the large-aggregation or latent correlation level: the realised correlation approached as \(\Delta\to\infty\), once clock-overlap losses and finite-response attenuation have disappeared. The coefficient \(\kappa_j\) is the effective relaxation rate with which the coupling-induced order flow in book \(j\) moves its reaction boundary in response to the directed spread \(z=p_1-p_2\). Since book 1 and book 2 move in opposite directions when the spread is nonzero, both motions reduce the same spread. Consequently the spread relaxation rate is the sum \(\kappa=\kappa_1+\kappa_2\).

4.3 Combined mechanism↩︎

When both effects are present, \[\rho_\Delta^{\mathrm{comb}} \simeq \rho_\infty \, \mathcal{A}_{12}(\Delta) \left[1-\frac{1-e^{-\kappa\Delta}}{\kappa\Delta}\right].\] For fractional clocks the exponential term is replaced by a Mittag–Leffler function [18]. See 10 for the details. The combined effect says that realised correlation builds up only when the two books both sample overlapping operational time and the coupling-induced reaction boundary response has had enough time to transmit the common movement.

5 Discussion↩︎

The contribution of this paper is that is combines three mature, but largely disjoint strands of the market microstructure literature: reaction–diffusion order book models, reduced–form explanations of the Epps effect, and discrete–event representations of high–frequency trading. Each of these strands is well developed in isolation.

Reaction–diffusion models successfully capture liquidity replenishment and price impact in single assets, but typically remain uncoupled across assets and make no claims about cross–correlations [2], [3], [19]. Conversely, Epps–effect models explain scale–dependent correlation using asynchronous sampling, lead–lag effects, or estimator bias, but do so at the level of observed prices rather than order–flow dynamics [9], [10], [17], [20]. The present work bridges these strands by embedding correlation emergence directly into a coupled reaction–diffusion representation of multiple order books.

The key idea is that correlation is not imposed exogenously through correlated noise or latent multivariate diffusions. Instead, it emerges endogenously from an explicit order–flow coupling between books. The coupling represents a pair–trader mechanism that responds to transient mispricing by placing side–selective orders near the reaction boundary. In contrast to Hawkes–process or multivariate diffusion models, where correlation parameters are fixed primitives [11], [12]; the effective correlation here is an output of the dynamics. This distinction is important: long–run correlation appears only after sufficient aggregation, while short–horizon correlation is suppressed because the coupling response is neither instantaneous nor perfectly synchronous.

A second important component is the role of the reaction boundary itself as the carrier of correlation. In reaction–diffusion order–book models, the transaction price is defined implicitly as a moving annihilation boundary between buy and sell densities. By analysing how the coupling perturbs the order density near this boundary, we can show that correlation transmission is filtered by the finite relaxation time of the reaction front. This provides a microstructural explanation for the Epps effect that we think is absent from earlier reduced–form derivations, where correlation attenuation is attributed solely to sampling effects or lead–lag structure [21], [22]. In the present framework, even perfectly synchronous observation does not eliminate short–horizon correlation decay if the coupling response itself is slow. This is important.

The paper also clarifies the relationship between discrete event dynamics and their continuum limits. Building on continuous–time random walk theory and fractional diffusion limits [6], [23], [24], we show how non–uniform event clocks and anomalous waiting times survive the diffusion limit and interact with coupling dynamics. This leads to a natural analytic decomposition of the Epps effect into three components: subordination–only attenuation due to asynchronous clocks, coupling–only attenuation due to finite response time, and a combined mechanism in which the coupling response is observed through non–uniform or fractional time changes. Earlier work separated asynchrony from lead–lag effects [10] and not endogenously bottom-up from the level of interacting order books.

Finally, the separation between the legacy discrete coupling used in numerical simulations and the regularised coupling employed in the analytic derivation is a minor contribution. The numerical model follows existing practice in agent–based [25] and order–book simulations [7], where piecewise rules govern order placement. For analytic tractability, we introduce a local regularisation that preserves the economic mechanism while removing singularities that obstruct linearisation around the symmetric state. This is important as this is structure that is strictly no pre-existing in the limiting processes required to move from the discrete event world to the continuum and should best be thought of as analogous to approximation functions required to impose continuity. Making this distinction explicit allows the numerical and analytic results to be compared without conflating discretisation artefacts with continuum behaviour. In this sense, the present work complements recent arguments that the Epps effect reflects the discrete event nature of markets [26] by showing precisely how such discreteness manifests after passing to the diffusion limit – and what additional assumptions are necessary.

6 Conclusion↩︎

By clarifying the discrete-to-continuum construction of coupled order books and using a regularised analytic representation of the pair-trader coupling, we obtained a transparent explanation of correlation emergence and the Epps effect. The analysis shows how asynchronous clocks, finite coupling response times, and their interaction each contribute distinct scale dependence.

7 Acknowledgements↩︎

We acknowledge conversations and some calculations being checked by Derick Diana, Byron Jacobs and Dominic Bauer. Copilot was used to help with some of the editing and proof reading, and to tidy up the references. All remaining errors are ours.

8 Subordination-only derivation↩︎

This appendix gives the derivation of the attenuation of realised correlation that arises solely from non-uniform observation or event clocks. Throughout this appendix the pair-trader coupling is switched off. The only mechanism reducing short-horizon correlation is imperfect overlap of the operational-time intervals sampled by the two assets.

8.1 Assumptions↩︎

We use the following assumptions.

1. There is a filtered probability space supporting two Brownian motions \((B^{(1)},B^{(2)})\) with instantaneous correlation \[d\langle B^{(1)},B^{(2)}\rangle_u=\rho_\infty\,du,\qquad -1\leq \rho_\infty\leq 1.\]

2. The latent operational-time prices are continuous semimartingales \[dX^{(j)}_u=\sigma_j\,dB^{(j)}_u,\qquad \sigma_j>0,\qquad j=1,2.\] Drifts may be added, but are omitted because their contribution to short-window covariance is lower order relative to the martingale covariance.

3. Calendar-time observations are generated by non-decreasing, right-continuous clocks \(T_j(t)\) independent of the Brownian motions. Their increments have finite first moments and, for the stationary formulae below, stationary increments.

4. Conditional on the clocks, the quadratic covariation of the two observed increments is obtained by the operational-time overlap of the two clock intervals.

The observed prices are \[p_j(t)=X^{(j)}_{T_j(t)},\qquad j=1,2.\] For a calendar interval \(I=(t,t+\Delta]\), define \[\Delta T_j(I)=T_j(t+\Delta)-T_j(t),\] and \[R^{(j)}_\Delta(t)=p_j(t+\Delta)-p_j(t).\]

8.2 Conditional covariance↩︎

Conditional on \(T_j\), the return of asset \(j\) is \[R^{(j)}_\Delta(t)=\sigma_j\int_{T_j(t)}^{T_j(t+\Delta)}dB^{(j)}_u.\] Hence, by the Itô isometry, \[\operatorname{Var}\!\left(R^{(j)}_\Delta(t)\mid T_j\right) =\sigma_j^2\int_{T_j(t)}^{T_j(t+\Delta)}du =\sigma_j^2\Delta T_j(I). \label{eq:sub-cond-var}\tag{11}\]

For the covariance, write the two stochastic integrals over the common operational-time axis using indicator functions: \[\begin{align} R^{(1)}_\Delta(t) &=\sigma_1\int_0^\infty \mathbf{1}_{(T_1(t),T_1(t+\Delta)]}(u)\,dB^{(1)}_u,\\ R^{(2)}_\Delta(t) &=\sigma_2\int_0^\infty \mathbf{1}_{(T_2(t),T_2(t+\Delta)]}(u)\,dB^{(2)}_u. \end{align}\] Therefore, \[\begin{align} &\operatorname{Cov}\!\left(R^{(1)}_\Delta(t),R^{(2)}_\Delta(t)\mid T_1,T_2\right)\nonumber\\ &\quad=\rho_\infty\sigma_1\sigma_2 \int_0^\infty \mathbf{1}_{(T_1(t),T_1(t+\Delta)]}(u) \mathbf{1}_{(T_2(t),T_2(t+\Delta)]}(u)\,du. \end{align}\] Define the sampled operational-time overlap \[\Theta_{12}(I) :=\left|(T_1(t),T_1(t+\Delta)]\cap (T_2(t),T_2(t+\Delta)]\right|,\] where \(|\cdot|\) denotes Lebesgue measure. Then \[\operatorname{Cov}\!\left(R^{(1)}_\Delta,R^{(2)}_\Delta\mid T_1,T_2\right) =\rho_\infty\sigma_1\sigma_2\Theta_{12}(I). \label{eq:sub-cond-cov}\tag{12}\] Taking expectations and using independence of the Brownian motions from the clocks gives \[\operatorname{Cov}\!\left(R^{(1)}_\Delta,R^{(2)}_\Delta\right) =\rho_\infty\sigma_1\sigma_2\mathbb{E}[\Theta_{12}(I)].\] Similarly, from 11 , \[\operatorname{Var}\!\left(R^{(j)}_\Delta\right) =\sigma_j^2\mathbb{E}[\Delta T_j(I)].\] Thus the realised correlation is \[\rho^{\rm sub}_\Delta =\rho_\infty\mathcal{A}_{12}(\Delta),\] for \[\mathcal{A}_{12}(\Delta) =\frac{\mathbb{E}[\Theta_{12}(I)]}{\sqrt{\mathbb{E}[\Delta T_1(I)]\mathbb{E}[\Delta T_2(I)]}}. \label{eq:sub-general-factor}\tag{13}\] This formula is exact under assumptions (A1)–(A4). All modelling approximations enter through the evaluation of the overlap factor \(\mathcal{A}_{12}\).

8.3 Poisson-refresh envelope↩︎

A tractable approximation is obtained by replacing the exact overlap process by a renewal envelope. Suppose the effective refresh processes of the two books are independent Poisson processes with rates \(\lambda_1\) and \(\lambda_2\). The probability that neither clock has refreshed over a lag \(u\) is \[\mathbb{P}(no refresh of either book over lag u)=e^{-(\lambda_1+\lambda_2)u}.\] Set \(\lambda_{12}=\lambda_1+\lambda_2\). Averaging the lag-dependent overlap over a window gives \[\begin{align} \frac{1}{\Delta}\int_0^\Delta\left(1-e^{-\lambda_{12}u}\right)du &=1-\frac{1}{\Delta}\int_0^\Delta e^{-\lambda_{12}u}du\\ &=1-\frac{1-e^{-\lambda_{12}\Delta}}{\lambda_{12}\Delta}. \end{align}\] Consequently, \[\rho^{\rm sub}_\Delta \simeq \rho_\infty \left[1-\frac{1-e^{-\lambda_{12}\Delta}}{\lambda_{12}\Delta}\right]. \label{eq:sub-poisson}\tag{14}\] The small-window expansion follows from \(e^{-x}=1-x+x^2/2+O(x^3)\): \[1-\frac{1-e^{-\lambda_{12}\Delta}}{\lambda_{12}\Delta} =\frac{\lambda_{12}\Delta}{2}+O(\Delta^2), \qquad \Delta\downarrow0.\] The large-window limit is \[1-\frac{1-e^{-\lambda_{12}\Delta}}{\lambda_{12}\Delta}\longrightarrow1, \qquad \Delta\to\infty.\] Thus asynchronous clocks suppress instantaneous correlation but recover the latent long-run correlation as aggregation increases.

8.4 Fractional-clock envelope↩︎

For heavy-tailed waiting times whose scaling limit is an inverse-stable clock of order \(0<\alpha\leq1\), exponential survival factors are replaced by Mittag–Leffler survival factors [18]. The fractional analogue of the previous envelope is \[\mathcal{A}_{12}^{(\alpha)}(\Delta) \simeq 1-\frac{1}{\Delta}\int_0^\Delta E_\alpha(-\lambda_{12}u^\alpha)\,du.\] Using the series definition \[E_\alpha(-\lambda u^\alpha)=\sum_{n=0}^{\infty}\frac{(-\lambda)^n u^{\alpha n}}{\Gamma(\alpha n+1)},\] and integrating term-by-term, \[\begin{align} \frac{1}{\Delta}\int_0^\Delta E_\alpha(-\lambda u^\alpha)\,du &=\sum_{n=0}^{\infty}\frac{(-\lambda)^n\Delta^{\alpha n}}{(\alpha n+1)\Gamma(\alpha n+1)}\\ &=\sum_{n=0}^{\infty}\frac{(-\lambda\Delta^\alpha)^n}{\Gamma(\alpha n+2)} =E_{\alpha,2}(-\lambda\Delta^\alpha). \end{align}\] Therefore \[\rho^{\rm sub,\alpha}_\Delta \simeq \rho_\infty\left[1-E_{\alpha,2}(-\lambda_{12}\Delta^\alpha)\right]. \label{eq:sub-fractional}\tag{15}\] When \(\alpha=1\), \(E_{1,2}(-x)=(1-e^{-x})/x\), and 15 reduces to 14 .

9 Coupling-only derivation↩︎

This appendix derives the Epps-type curve generated by the finite response time of the pair-trader coupling when calendar-time clocks are synchronous. The purpose is to show explicitly how a local order-flow coupling near the moving reaction boundary produces a mean-reverting price-level approximation and hence a scale-dependent realised correlation.

9.1 Assumptions↩︎

We assume the following.

1. For each book \(j\), the signed order density \(\varphi^{(j)}(x,t)\) is continuously differentiable in \(t\) and twice continuously differentiable in \(x\) in a neighbourhood of the reaction boundary.

2. The transaction price \(p_j(t)\) is the unique transverse zero crossing of \(\varphi^{(j)}\): \[\varphi^{(j)}(p_j(t),t)=0, \qquad \mathcal{L}_j(t):=\partial_x\varphi^{(j)}(p_j(t),t)\neq0.\]

3. The pair-trader source is local near the reaction boundary and has finite first spatial moment.

4. The coupling is weak enough that the slope \(\mathcal{L}_j(t)\) and the unperturbed book shape may be frozen over the short response calculation.

5. The idiosyncratic price innovations generated by uncoupled order-flow shocks are represented, after the local reduction, by Brownian terms with volatilities \(\sigma_j>0\).

9.2 Zero-crossing reduction↩︎

The moving boundary is defined implicitly by \[F_j(p_j(t),t):=\varphi^{(j)}(p_j(t),t)=0.\] Differentiating gives \[0=\frac{d}{dt}\varphi^{(j)}(p_j(t),t) =\partial_t\varphi^{(j)}(p_j(t),t) +\partial_x\varphi^{(j)}(p_j(t),t)\dot{p}_j(t).\] Hence \[\dot{p}_j(t) =-\frac{\partial_t\varphi^{(j)}(p_j(t),t)}{\mathcal{L}_j(t)}. \label{eq:zc-price-velocity}\tag{16}\] If the order-density equation contains a source term \(c^{(j,k)}=s^{(j)}+\ell^{(j,k)}+\delta^{(j)}\), then the immediate contribution of the coupling to 16 is formally \[\dot{p}_{j,\rm coup}(t) =-\frac{\ell^{(j,k)}(p_j(t),t)}{\mathcal{L}_j(t)}.\] However, the analytic coupling used here is odd at the mid-price and satisfies \(q_j(0)=0\). Therefore its point value at the zero crossing vanishes. The price displacement is consequently captured by the first non-zero spatial moment of the source around the crossing, equivalently by projecting the source perturbation onto the translational mode of the reaction front. In the frozen-slope approximation this gives \[\dot{p}_{j,\rm coup}(t) \simeq -\frac{1}{\mathcal{L}_j(t)} \int_{\mathbb{R}} y\,\ell^{(j,k)}(p_j(t)+y,t)\,dy, \label{eq:moment-projection}\tag{17}\] up to a convention-dependent positive normalising constant that may be absorbed into the definition of the effective coupling strength. The frozen-slope approximation is a local sharp-interface approximation: over the short response calculation the reaction front is displaced, but its local shape and slope \(L_j(t)=\partial_x\phi^{(j)}(p_j(t),t)\) are treated as fixed. This is the same type of local projection used in sharp-interface and level-set reductions of diffuse field models to moving-boundary laws [27].

9.3 Regularised analytic coupling↩︎

Let \[y_j=x-p_j(t), \qquad z_{jk}=p_j(t)-p_k(t).\] The regularised local-response coupling is \[\ell^{(j,k)}(x,t) =\gamma_{jk}z_{jk}(t)q_j(y_j) W(y_j,z_{jk};\varepsilon), \label{eq:coup-reg-local}\tag{18}\] with \(\gamma_{jk}\geq0\) and \[q_j(y)=-\lambda_j\mu_j y e^{-\mu_j y^2}, \qquad \lambda_j>0,\quad \mu_j>0.\] A smooth side-selection is \[W(y,z;\varepsilon)=\frac{1}{2}\left[1+\tanh\!\left(\frac{yz}{\varepsilon}\right)\right], \qquad \varepsilon>0.\] In the hard-side limit one obtains \[W_0(y,z)=\mathbf{1}_{\{yz>0\}}.\] The crucial point is that 18 is finite and differentiable as \(z_{jk}\to0\). This property is not shared by the legacy piecewise coupling containing inverse powers of \(z_{jk}\). The regularisation is therefore not merely notational: it is required for the Taylor expansion around the symmetric state \(p_1=p_2\).

9.4 Evaluation of the first moment↩︎

For \(z_{jk}>0\), the hard-side limit selects \(y>0\). Define \[M_j^{+}:=\int_{0}^{\infty}yq_j(y)\,dy.\] With \(q_j(y)=-\lambda_j\mu_j y e^{-\mu_j y^2}\), \[\begin{align} M_j^{+} &=-\lambda_j\mu_j\int_0^\infty y^2e^{-\mu_j y^2}\,dy. \end{align}\] Using \[\int_0^\infty y^2e^{-\mu y^2}\,dy =\frac{\sqrt\pi}{4\mu^{3/2}},\] we obtain \[M_j^{+}=-\frac{\lambda_j\sqrt\pi}{4\sqrt{\mu_j}}<0. \label{eq:moment-explicit}\tag{19}\] For \(z_{jk}<0\), the selected side is \(y<0\). Since \(yq_j(y)\) is even, the same numerical moment is obtained. Thus the moment contribution is proportional to \(z_{jk}\) and has the sign required for mean reversion after choosing the orientation of the signed order density consistently.

Substituting 18 into 17 gives \[\dot{p}_{j,\rm coup}(t) \simeq -\frac{\gamma_{jk}z_{jk}(t)}{\mathcal{L}_j(t)}M_j,\] where \(M_j\) denotes the side-selected first moment with the sign convention of book \(j\). We define the positive effective relaxation coefficient by \[\kappa_j:=\frac{\gamma_{jk}|M_j|}{|\mathcal{L}_j|}>0,\] so that, for the two-book spread \(z=p_1-p_2\), \[\begin{align} dp_1(t)&=\sigma_1dW_1(t)-\kappa_1z(t)dt,\tag{20} \\ dp_2(t)&=\sigma_2dW_2(t)+\kappa_2z(t)dt. \tag{21} \end{align}\] This is the local price-level approximation induced by the regularised coupling.

9.5 Spread dynamics↩︎

Subtracting the two equations 20 and 21 gives \[dz(t)=-\kappa z(t)dt+\sigma_1dW_1(t)-\sigma_2dW_2(t), \qquad \kappa:=\kappa_1+\kappa_2.\] If \(d\langle W_1,W_2\rangle_t=\rho_0dt\), then the spread-noise variance is \[\sigma_z^2=\sigma_1^2+\sigma_2^2-2\rho_0\sigma_1\sigma_2.\] Thus \[dz(t)=-\kappa z(t)dt+\sigma_zdW_z(t),\] which is an Ornstein–Uhlenbeck process with response time \(\kappa^{-1}\).

9.6 Covariance build-up↩︎

The finite response time is represented by an even cross-covariance kernel \[c_{12}(u)=C\frac{\kappa}{2}e^{-\kappa|u|}, \label{eq:exp-cov-kernel}\tag{22}\] where \(C\) fixes the long-window covariance scale. The covariance of returns over a window of length \(\Delta\) is the convolution of this kernel with the triangular overlap weight: \[\begin{align} \operatorname{Cov}_\Delta &=\int_{-\Delta}^{\Delta}(\Delta-|u|)c_{12}(u)\,du\\ &=C\kappa\int_0^\Delta(\Delta-u)e^{-\kappa u}\,du. \end{align}\] Now \[\begin{align} \int_0^\Delta(\Delta-u)e^{-\kappa u}\,du &=\Delta\int_0^\Delta e^{-\kappa u}\,du -\int_0^\Delta ue^{-\kappa u}\,du\nonumber\\ &=\Delta\frac{1-e^{-\kappa\Delta}}{\kappa} -\left[-\frac{\Delta e^{-\kappa\Delta}}{\kappa} +\frac{1-e^{-\kappa\Delta}}{\kappa^2}\right]\nonumber\\ &=\frac{\Delta}{\kappa}-\frac{1-e^{-\kappa\Delta}}{\kappa^2}. \end{align}\] Therefore \[\operatorname{Cov}_\Delta =C\left[\Delta-\frac{1-e^{-\kappa\Delta}}{\kappa}\right].\] If the variances satisfy \[\operatorname{Var}(R^{(j)}_\Delta)\simeq v_j\Delta, \qquad v_j>0,\] then \[\rho^{\rm coup}_\Delta \simeq \rho_\infty \left[1-\frac{1-e^{-\kappa\Delta}}{\kappa\Delta}\right], \qquad \rho_\infty:=\frac{C}{\sqrt{v_1v_2}}. \label{eq:coup-final}\tag{23}\] At small scales, \[\rho^{\rm coup}_\Delta \sim \rho_\infty\frac{\kappa\Delta}{2}, \qquad \Delta\downarrow0,\] and at large scales, \[\rho^{\rm coup}_\Delta\to \rho_\infty, \qquad \Delta\to\infty.\]

10 Combined derivation↩︎

This appendix combines the two attenuation mechanisms derived separately above. The first mechanism is incomplete overlap of event or observation clocks. The second is finite relaxation of the coupling-induced response at the moving reaction boundary. The combined formulae below are approximations unless the clock process and the coupling response are independent and enter multiplicatively.

10.1 Assumptions↩︎

We assume the following.

1. In operational time, the coupling-induced covariance response has a completely monotone relaxation survival function \(h(u)\) with response density \(K(u)=-h'(u)\).

2. The observation clocks are independent of the idiosyncratic price innovations and independent of the local coupling kernel to first order.

3. Variance growth remains approximately linear over the aggregation scales considered, or else the final correlation should be interpreted with the corresponding variance normalisation explicitly included.

4. The product formula is a first-order separable approximation. It is not exact when the clock process changes the local coupling intensity or when pair-trader activity is itself triggered by the same clock events that generate observations.

10.2 Separable and light-tailed↩︎

When clocks attenuate the amount of common operational time sampled, while the coupling response is approximately exponential in calendar time, the covariance scale factor separates as \[\rho_\Delta^{\rm comb} \simeq \rho_\infty\mathcal{A}_{12}(\Delta)F_\kappa(\Delta), \label{eq:comb-separable}\tag{24}\] where \[F_\kappa(\Delta) =1-\frac{1-e^{-\kappa\Delta}}{\kappa\Delta}.\] Thus, for the Poisson-refresh envelope, \[\rho_\Delta^{\rm comb} \simeq \rho_\infty \left[1-\frac{1-e^{-\lambda_{12}\Delta}}{\lambda_{12}\Delta}\right] \left[1-\frac{1-e^{-\kappa\Delta}}{\kappa\Delta}\right]. \label{eq:comb-poisson}\tag{25}\] Equation 25 should be read as a leading-order factorisation. It is most appropriate when the clock process determines which increments are compared, while the local response of the reaction boundary is governed by an independent relaxation timescale \(\kappa^{-1}\).

10.3 Fractional relaxation↩︎

If the coupling relaxation is subordinated by an inverse-stable clock of index \(0<\alpha\leq1\), exponential relaxation is replaced by Mittag–Leffler relaxation. In operational time the survival function solves [18] \[\frac{dh}{du}=-\kappa h(u), \qquad h(0)=1,\] with \(h(u)=e^{-\kappa u}\). Under inverse-stable subordination, the calendar-time survival function \(h_\alpha(t)\) solves the Caputo fractional relaxation equation \[{}^C D_t^\alpha h_\alpha(t)=-\kappa h_\alpha(t), \qquad h_\alpha(0)=1.\] Taking Laplace transforms gives \[s^\alpha\widehat h_\alpha(s)-s^{\alpha-1} =-\kappa\widehat h_\alpha(s),\] so that \[\widehat h_\alpha(s)=\frac{s^{\alpha-1}}{s^\alpha+\kappa}.\] Inverting yields \[h_\alpha(t)=E_\alpha(-\kappa t^\alpha).\] The corresponding response density is \[\begin{align} K_{\alpha,\kappa}(t) &=-\frac{d}{dt}E_\alpha(-\kappa t^\alpha)\\ &=\kappa t^{\alpha-1}E_{\alpha,\alpha}(-\kappa t^\alpha), \qquad t>0. \end{align}\] This follows by differentiating the series for \(E_\alpha\) term by term.

10.4 Fractional covariance build-up↩︎

The normalised build-up factor over an aggregation window is \[F_\alpha(\Delta) =\frac{1}{\Delta}\int_0^\Delta(\Delta-u)K_{\alpha,\kappa}(u)\,du.\] Since \(K=-h'\), integration by parts gives \[\begin{align} \int_0^\Delta(\Delta-u)K(u)\,du &=-\int_0^\Delta(\Delta-u)h'(u)\,du\\ &=\Delta h(0)-\int_0^\Delta h(u)\,du, \end{align}\] where the endpoint term at \(u=\Delta\) vanishes. Because \(h(0)=1\), \[F_\alpha(\Delta) =1-\frac{1}{\Delta}\int_0^\Delta E_\alpha(-\kappa u^\alpha)\,du.\] Using the identity \[\frac{1}{\Delta}\int_0^\Delta E_\alpha(-\kappa u^\alpha)\,du =E_{\alpha,2}(-\kappa\Delta^\alpha),\] we obtain \[F_\alpha(\Delta)=1-E_{\alpha,2}(-\kappa\Delta^\alpha). \label{eq:fractional-build-factor}\tag{26}\] Therefore the fractional coupling-only realised correlation is \[\rho^{(\alpha)}_\Delta \simeq \rho_\infty\left[1-E_{\alpha,2}(-\kappa\Delta^\alpha) \right]. \label{eq:fractional-correlation}\tag{27}\] For \(\alpha=1\), \[E_{1,2}(-x)=\frac{1-e^{-x}}{x},\] so 27 reduces to the ordinary coupling curve.

10.5 Small-scale asymptotics↩︎

The two-parameter Mittag–Leffler expansion [18] is \[E_{\alpha,2}(-\kappa\Delta^\alpha) =\frac{1}{\Gamma(2)}- \frac{\kappa\Delta^\alpha}{\Gamma(\alpha+2)}+O(\Delta^{2\alpha}).\] Since \(\Gamma(2)=1\), \[F_\alpha(\Delta) =\frac{\kappa\Delta^\alpha}{\Gamma(\alpha+2)}+O(\Delta^{2\alpha}).\] Thus fractional time changes the initial build-up from order \(\Delta\) to order \(\Delta^\alpha\). For \(0<\alpha<1\), the build-up is slower in calendar time at short horizons, reflecting heavy-tailed waiting times.

10.6 Combined approximation↩︎

If both clock-overlap attenuation and fractional relaxation are retained, the first-order separable formula becomes \[\rho^{\rm comb,\alpha}_\Delta \simeq \rho_\infty\mathcal{A}_{12}^{(\alpha_c)}(\Delta) \left[1-E_{\alpha_r,2}(-\kappa\Delta^{\alpha_r})\right],\] where \(\alpha_c\) describes the observation-clock overlap and \(\alpha_r\) describes the relaxation response. In the simplest single-clock scaling one may set \(\alpha_c=\alpha_r=\alpha\), but this equality is a modelling assumption rather than a mathematical necessity.

11 Coupling Equivalence↩︎

We refer the interested reader to [14], [28]–[30]. Briefly, the distinction between pointwise and weak equivalence seems accepted or at least standard in regularised-interface and immersed-boundary calculations. A singular or discontinuous object is first replaced by a bounded smooth kernel, and the comparison with the sharp object is then made only after testing against the modes or moments that enter the reduced equation. Our use of \(W(y,z;\varepsilon)\) should be read in this sense. The regularised pair-trader source is not required to coincide pointwise with the legacy piecewise source near \(z=0\), where the latter contains singular \(1/z\)-type branches. What is required for the price-level reduction is only that the projection of the source onto the reaction-front displacement mode. This seems analogous to sharp-interface asymptotics for diffuse-interface models, where a finite-width interface is introduced first and the reduced moving-boundary law is obtained by projection before taking the sharp-interface limit [28].

Here the two couplings are not pointwise identical. The idea here is that they are equivalent only in a weak, local, first-moment sense after imposing a regularisation and taking a controlled small-spread, small-grid limit.

11.1 Objects being compared↩︎

Let \(p_j\) and \(p_k\) be the reaction-boundary prices of two books, and write \[y=x-p_j, \qquad z=z_{jk}=p_j-p_k.\] The numerical coupling is piecewise and side-selective. In schematic form it applies a source shape \(g_j\) to one side of the book, with branches depending on the signs of \(y\) and \(z\). Some legacy branches contain inverse powers of \(z\), or scaled arguments such as \(y/z\). These terms are harmless away from \(z=0\) on a finite grid, provided the simulation never evaluates the singular branch at exact equality, but they obstruct analytic expansion about the symmetric state \(z=0\).

The analytic coupling is \[\ell_{\rm reg}^{(j,k)}(x,t) =\gamma_{jk}zq_j(y)W(y,z;\varepsilon), \label{eq:eq-reg}\tag{28}\] where \(q_j\) is odd and locally linear at \(y=0\), and \(W\) is bounded with hard-side limit \[W_0(y,z)=\mathbf{1}_{\{yz>0\}}.\]

11.2 Why pointwise equivalence fails↩︎

Pointwise equivalence near \(z=0\) would require \[\ell_{\rm disc}(y,z)\to \ell_{\rm reg}(y,z)\] for fixed \(y\) as \(z\to0\). This cannot hold for the legacy expression if any branch contains \(y/z\) or \(1/z\), because such terms do not have a finite pointwise limit at \(z=0\). For example, if a branch contains \(g_j(y/z)\) and \[g_j(\eta)=-\lambda_j\mu_j\eta e^{-\mu_j\eta^2},\] then for fixed \(y\neq0\), \[g_j(y/z)=-\lambda_j\mu_j\frac{y}{z}\exp\!\left[-\mu_j\frac{y^2}{z^2}\right] \longrightarrow0\] for Gaussian decay, while with the finite-domain convention \(\exp(+\mu_j y^2/z^2)\) the expression diverges. In either case the limit is not the linear local response \(zq_j(y)\). Hence regularisation is mathematically necessary before taking a Taylor expansion in \(z\).

11.3 Weak local equivalence criterion↩︎

The price equation depends on the coupling through its projection onto the reaction-front displacement mode. In the frozen-slope approximation this projection is represented by the first moment \[\mathcal{M}[\ell](z)=\int_{\mathbb{R}}y\ell(y,z)\,dy.\] We therefore define weak local equivalence as follows. The discrete and regularised couplings are equivalent at the level relevant for the price reduction if \[\mathcal{M}[\ell_{\rm disc}^{\Delta x,\varepsilon}](z) =\eta_j z+o(z)\] and \[\mathcal{M}[\ell_{\rm reg}^{\varepsilon}](z) =\eta_j z+o(z) \label{eq:weak-equiv}\tag{29}\] with the same finite coefficient \(\beta_j\), after matching the discrete coupling strength to \(\gamma_{jk}\).

For the regularised coupling, \[\begin{align} \mathcal{M}[\ell_{\rm reg}](z) &=\gamma_{jk}z\int_{\mathbb{R}}yq_j(y)W(y,z;\varepsilon)\,dy. \end{align}\] In the hard-side limit, for \(z>0\), \[\mathcal{M}[\ell_{\rm reg}](z) =\gamma_{jk}z\int_0^\infty yq_j(y)\,dy =\gamma_{jk}M_j^+z.\] Thus the regularised coupling satisfies 29 with \[\beta_j=\gamma_{jk}M_j^+.\] For the Gaussian source shape, \[M_j^+=-\frac{\lambda_j\sqrt\pi}{4\sqrt{\mu_j}}.\] The discrete coupling has the same continuum price-level effect if its lattice first moment satisfies \[\Delta x\sum_i (x_i-p_j)\ell_{{\rm disc},i}^{(j,k)} =\gamma_{jk}M_j^+z+o(z)+o(1) \qquad \label{eq:lattice-moment-match}\tag{30}\] where as \(\Delta x\to0\), then \(z\to0\), or the analogous joint limit with a specified relation between \(z\), \(\Delta x\), and \(\varepsilon\).

11.4 The regularisation parameter↩︎

The smooth selector \[W(y,z;\varepsilon)=\frac{1}{2}\left[1+\tanh\left(\frac{yz}{\varepsilon}\right)\right]\] introduces a transition layer of width approximately \(\varepsilon/|z|\) in the \(y\) variable. If \(\varepsilon\) is fixed and \(z\to0\), then \(W\to1/2\) pointwise, and because \(yq_j(y)\) is even, \[\int_{\mathbb{R}}yq_j(y)W(y,z;\varepsilon)dy \to \frac{1}{2}\int_{\mathbb{R}}yq_j(y)dy =\int_0^\infty yq_j(y)dy.\] Thus the first-moment coefficient remains finite and equals the hard-side coefficient for symmetric kernels. If an asymmetric kernel is used, the order of limits \(\varepsilon\to0\) and \(z\to0\) must be specified explicitly.

A practical matching condition is \[\Delta x\ll \sqrt{\varepsilon}\ll L_q,\] where \(L_q\) is the spatial scale on which \(q_j\) varies. The first inequality makes the smoothed side selector resolvable on the lattice, while the second keeps the selector local relative to the order-book source shape.

11.5 Impact on the diffusion limit↩︎

The uncoupled lattice update converges to a reaction–diffusion or fractional reaction–diffusion equation under the scaling \[D_\alpha=\lim_{\Delta x\to0,\,\Delta t\to0} \frac{r}{2}\frac{\Delta x^2}{\Delta t^\alpha}.\] For the coupled equation to have a finite continuum limit, the source term must remain bounded in the topology used for convergence. The regularised coupling satisfies this requirement because \[|\ell_{\rm reg}(y,z)| \leq \gamma_{jk}|z|\,|q_j(y)|,\] and \(q_j\) is integrable with finite first moment. Therefore the coupling enters the continuum equation as an ordinary source term of order \(O(z)\).

By contrast, an unregularised branch containing \(1/z\) or \(y/z\) can fail to be uniformly bounded as \(z\to0\). If the spread \(z\) is itself of order \(\Delta x\) near the continuum limit, such a branch can scale like \(O(1/\Delta x)\) or worse. This can dominate the diffusive scaling and produce a singular source measure rather than the intended finite reaction term. Hence the diffusion limit and the symmetric-spread limit do not commute for the unregularised coupling.

The appropriate limiting order is therefore:

1. introduce a bounded regularised coupling at the lattice level;

2. take the diffusion limit \(\Delta x,\Delta t\to0\) with \(D_\alpha\) fixed and with source moments uniformly bounded;

3. perform the local small-spread expansion of the continuum source to obtain the finite response coefficient \(\kappa_j\);

4. if desired, take a hard-side limit of the selector after the first moment has been matched.

This equivalence is sufficient for the derivation of the Epps-effect curve, because that derivation depends only on the finite first-order response rate \(\kappa=\kappa_1+\kappa_2\), not on pointwise equality of the full source fields.

References↩︎

---

1. The legacy implementation contains branches involving inverse powers of \(\Delta p_{jk}\) and the expression \(j/\Delta p_{jk}\), where \(j\) is also the order-book index. The regularised notation used here is not intended to change the economic mechanism of the simulations; it is a local analytic replacement that is well-defined as \(p^{(j)}\to p^{(k)}\). For exact numerical reproducibility of [8] one should use the piecewise coupling 4 , while the analytic reduction should use the regularised coupling in 7 .↩︎

2. If the finite-domain numerical convention \(\exp(+\mu_jy^2)\) is retained, our regularised formulation remains unchanged after replacing \(q_j\) by \(-\lambda_j\mu_jy\exp(+\mu_jy^2)\).↩︎