\section{Introduction}

In early 2012 the journal \textit{Nature} had a News Feature article on ``tough science'' listing five experiments ``as hard as finding the Higgs'' \cite{Jones_Nature2012}. One experiment in the list attempts to provide the first observation of parity violation (PV), associated with the weak force, in molecular systems using high-resolution laser spectroscopy\cite{Crassous_OBC2005}.
It is an example of a spectroscopic test of fundamental physics in the low-energy regime, using atoms and molecules, in vivid contrast to the Large Hadron Collider, 27 km long and operating at several TeV \cite{Erler_PPNP2005,Safronova_RevModPhys.90.025008}.

It was Lee and Yang who in 1956 pointed out that there was insufficient experimental evidence for the conservation of parity in processes involving the
weak interaction and proposed experiments to test possible parity violation\cite{Lee_Yang_PhysRev.104.254}, soon to be confirmed by experiments\cite{wu1957experimental,Postma_Physica1957,Garwin_PhysRev.105.1415,Friedman_PhysRev.105.1681.2}. The processes that were studied are mediated by \textit{charged} vector bosons, in contrast to the neutral photon of electromagnetic interactions, and are therefore incompatible with stable atoms and molecules.
The situation changed with the development of electroweak theory\cite{Glashow_RevModPhys.52.539,Weinberg_RevModPhys.52.515,Salam_RevModPhys.52.525} and in particular the prediction of the existence of a neutral partner, $Z^0$, of the charged $W^\pm$ bosons. Early studies about the feasibility of observing parity violation in stable atoms were rather pessimistic\cite{Zeldovich_JETP1959,Michel_PR1965}, but limited attention to the hydrogen atom. M. A. and C. Bouchiat noted that PV has a strong scaling with nuclear charge and therefore suggested to look for such effects in heavy atoms\cite{Bouchiat_PLB1974,Bouchiat_JPF1974} and, indeed, the first observation of atomic PV was in the form of optical activity in a vapor of bismuth atoms\cite{Barkov_JETPlett1978,Barkov_PLB1979}. Since then PV has been observed in several heavy atoms:  thallium\cite{Conti_PRL1979,Bucksbaum_PRL1981,Drell_PRL1984,Vetter_PRL1995,Edwards_PRL1995},  cesium\cite{Bouchiat_PLB1982,Bouchiat_JPF1986,wood1997measurement,Bennett_PRL1999,Guena_PRL2003,Guena_PRA2005},  lead\cite{Emmons_PRL1983,Meekhof_PRL1993,Phipp_JPB1996}, bismuth\cite{Macpherson_PRL1991} and ytterbium\cite{Tsigutkin_PRL2009,Antypas_NaturePhysics2019}. For a recent review on atomic parity violation, see Ref.~\citenum{Roberts_ARNPS2015}.

For chiral molecules parity violation induces a finite, albeit minute energy difference between enantiomers\cite{letokhov1975difference,zel1977energy,Rein_PLA1979}, making them formally diastereomers, an intriguing fact in view of the observed chiral preferences in biological systems\cite{Bonner_TSC1988,Quack_Ang2002,Tsarev2009PPN}, but yet to be observed. The experiment under construction in Paris will search for the signature of parity violation in the form of a frequency shift in the vibrational spectra of left- and right-handed molecules. It aims for a measurement precision below 0.1 Hz and is based on the method of the Doppler-free two-photon Ramsey fringes in a supersonic molecular beam\cite{Crassous_OBC2005,Darquie:2010:chirality}. Recent progress is reported in Ref.~\citenum{Cournol_QE2019}, whereas information about other experiments and  molecular (as well as atomic) PV in general is found in an excellent review by Berger and Stohner\cite{Berger2019WCM}.

The role of theory at this point is to guide experiment. We have developed a computational protocol for the rapid screening of candidate molecules, based on 2-component relativistic density functional theory (DFT)\cite{DeMontigny_PCCP2010}. However, calculations of higher accuracy are needed, in the short run to calibrate the current computational protocol for molecular PV calculations, in the longer run to confront a possibly successful experiment. In this context coupled cluster theory stands out as the method of choice\cite{Bartlett_Stanton_RCC1994,Crawford_Schaefer_RCC2000,Bartlett_RevModPhys.79.291,Liu_Cheng_WIRE2021}. In 2000 Thyssen \textit{et al.} reported 4-component wavefunction-based correlated calculations of the parity-violating energy of the \ce{H2O2} and \ce{H2S2} molecules using a finite-field
approach at the level of 2nd-order M{\o}ller--Plesset perturbation theory (MP2), coupled-cluster singles-and-doubles (CCSD), CCSD with approximate triples correction (CCSD(T)) as well as configuration-interaction singles-and doubles (CISD)\cite{Thyssen_PRL2000}. The \ce{H2X2}-series of molecules are not good candidates for a PV experiment due to the almost free rotation around the X--X bond, but have been much used for benchmarking and analysis purposes due to their simple structures (see for instance Refs.\cite{mason1984parity,Bakasov_JCP_1998,lazzeretti1997calculation,lazzeretti1997calculation,laerdahl1999fully,Hennum_CPL2002,Berger_vanWullen_JCP2005,Berger_JCP2008,bast2011analysis,Horny_MP2015}). In 2005 van Stralen \textit{et al.} reported the first implementation of analytic one-electron
properties at the 4-component relativistic MP2 level and reported PV energies for \ce{H2X2}, (X=O,S,Se,Te) \cite{van2005first}. More recently, in 2016, these analytic calculations were extended to the CCSD level by Shee \textit{et al.} \cite{shee2016analytic}.

In the present work we repeat the calculations of Shee \textit{et al.} \cite{shee2016analytic}, but with a number of improvements: i) The coupled cluster calculations have been carried out using a new coupled cluster code, ExaCorr, implemented for massively parallel GPU-accelerated computing architectures\cite{pototschnig2021implementation}. ii) We have extended the Dyall-type all-electron basis sets with tight exponents optimized for PV calculations. iii) Using \ce{H2Se2}, as detailed in the Supplemental Material, we have investigated the effect of truncation of the occupied and virtual spaces  on the PV energy. iv) We have extended the calculations to include \ce{H2Po2}. v) Most importantly, we have to our knowledge carried out the first study of the effect of quantum electrodynamics (QED) on molecular PV energies.

We are here interested in the electron self-energy (SE) and the vacuum polarization (VP), both at the origin of the Lamb shift\cite{lamb1947fine}. QED, as expressed through the scattering matrix formalism, can describe atoms and molecules to amazing precision, but is in practice limited to few-electron systems, in particular due to the slow convergence of electron correlation within this formalism\cite{Indelicato_JPhysB2019,mohr:QEDrev,Sapirstein_PhysRevA2002,Sapirstein_PhysRevA.91.062508,Korobov_PhysRevLett2016,Korobov_PhysRevLett2017,Yerokhin_PhysLettA1995}. In recent years there has been an increased interest in the use of effective QED potentials\cite{shepler2005ab,Thierfelder2010PRA,Artemyev_chap2016,pavsteka2017relativistic}. An early example of such a potential is the Uehling potential for vacuum polarization\cite{Uehling1935}. The electron self-energy is more complicated to represent by a potential due to its delocal nature, but a number of such potentials are now available\cite{Pyykko2003a,Flambaum2005,Shabaev_PhysRevA.88.012513}. QED-effects have also been included in the fitting of relativistic effective core potentials\cite{hangele2012accurate,hangele2013relativistic}. We have implemented effective QED potentials in the DIRAC code for general relativistic molecular calculations\cite{saue2020dirac}. A full account of our implementation will be given elsewhere\cite{Sunaga2021qed}; very recently, a similar implementation has been reported by Skripnikov\cite{Skripnikov_JCP2021}. QED-effects on parity violation have previously been studied in  atoms and then with focus on transition amplitudes\cite{Blundell_PRD1992,johnson2001vacuum,milstein2002radiative,Shabaev_PRA2005,kuchiev2003radiative,Roberts_PRA2013a,Roberts_PRA2013b,Sahoo_MP2017,Yu_Khetselius_2017}. We believe that the present work is the first to study radiative corrections to molecular PV energies.

Our paper is organized as follows: In Section \ref{seq:effQED} we explain how effective QED potentials are included in our calculations. Section \ref{sec:PV} gives a brief introduction to
parity violation in chiral molecules. Computational details are given in in Section \ref{sec:comp}. In Section \ref{sec:results} we present and discuss our results, before concluding in Section \ref{sec:conc}. We employ SI-based atomic units throughout this paper.
\section{Theory}
\subsection{Relativistic Hamiltonian with effective QED potentials}\label{seq:effQED}

Our calculations are based on the Dirac--Coulomb Hamiltonian
\begin{equation}
{\hat H}_{\rm{DC}} = V_{nn} + \sum\limits_i {\hat h}(i) + \frac{1}{2}\sum\limits_{i \ne j}\frac{1}{r_{ij}}
\end{equation}
Here, the first and third terms represent the nucleus-nucleus and electron-electron repulsion, respectively. The one-electron part is given by the Dirac Hamiltonian
\begin{equation}
{\hat h} =   {\beta}m_{\rm{e}}c^2 + c\left(\boldsymbol{\alpha}\cdot{\bf p}\right) + V_{\rm{en}};\quad V_{\rm{en}}=-e\phi_n
\end{equation}
in the scalar potential $\phi_n$ of fixed nuclei and will be extended by effective QED-potentials $V_{\rm{effQED}}$
\begin{equation}
{\hat h} \rightarrow {\hat h} + V_{\rm{effQED}};\quad V_{{\rm{effQED}}} = \sum\limits_A\left(V^{\rm{VP}}_A + V^{\rm{SE}}_A\right)
\end{equation}
where $V^{\rm{VP}}_A$ and $V^{\rm{SE}}_A$ refer to the contributions from the vacuum polarization and the electron self-energy, respectively, associated with nucleus $A$.
We employ the Uehling potential \cite{Uehling1935} for the vacuum polarization term $V^{\rm{VP}}$.  The Uehling potential for a spherically symmetric extended nucleus can be expressed as follows \cite{Fullerton1976}
\begin{equation}
\label{eqn:Uehling_exact}
V^{{\rm{VP}}}_A(r) =  -\frac{2\alpha{\mkern0.75mu\mathchar '26\mkern -9.75mu\lambda}_e}{3}\frac{Z_Ae^2}{r}\int^{\infty}_0dr'r'\rho_A\left(r'\right)\left[K_0\left(\frac{2}{{\mkern0.75mu\mathchar '26\mkern -9.75mu\lambda}_e}\left|r-r'\right|\right)\right.  \left.-K_0\left(\frac{2}{{\mkern0.75mu\mathchar '26\mkern -9.75mu\lambda}_e}\left|r+r'\right|\right)\right].
\end{equation}
where
\begin{equation}
K_0(x)=\int^{\infty}_1 dt\exp{[-xt]}\left(\frac{1}{t^3}+\frac{1}{2t^5}\right)\left(t^2-1\right)^{1/2},
\end{equation}
$\alpha$ the fine structure constant, ${{\mkern0.75mu\mathchar '26\mkern -9.75mu\lambda}}_e$ the reduced Compton wavelength and $\rho_A$ is the normalized charge density of nucleus $A$.

For the self-energy term $V_{{\rm{SE}}}$, we have employed two kinds of effective potentials. The first one is the Gaussian model potential proposed by Pyykk\"{o} and Zhao (PZ) \cite{Pyykko2003a},
\begin{equation}
{V^{{\rm{SE}}}_A} \rightarrow {V^{\rm{PZ}}_A}\left( r \right) = B(Z_A)\exp \left( { - \beta (Z_A){r^2}} \right),
\end{equation}
where the parameters $B$ and $\beta$ are second-order polynomials in nuclear charge with coefficients obtained by fitting to SE data ($2s_{1/2}$ SE energy shift of hydrogen-like ions\cite{Indelicato_PhysRevA.58.165,Beier_PRA1998} and relative SE shift of hyperfine transition energies of lithium-like ions\cite{Boucard_EPJD2000}) in the range $29\le Z \le 83$.

The second effective SE potential was proposed by Flambaum and Ginges (FG) \cite{Flambaum2005}. It is formulated in the language of form factors expressing the effective interaction of an electron with the electromagnetic field, radiative corrections included\cite{Schweber_RQFT,Berestetskii:QED1982}. The potential is written as a sum of three contributions
\begin{equation}
{V^{{\rm{SE}}}_A} \rightarrow {V^{\rm{FG}}_A}\left( r \right)= V^\textrm{MAG}_A + V^\textrm{HF}_A + V^\textrm{LF}_A,
\end{equation}
where the first term is associated with the magnetic form factor
\begin{equation}\label{eqn:MAG}
V^\textrm{MAG}_A(\mathbf{r})= \frac{\alpha{\mkern0.75mu\mathchar '26\mkern -9.75mu\lambda}_e}{4\pi} i \beta\boldsymbol{\alpha} \cdot \nabla \left[\phi_{n;A}(r) \left(\int^{\infty}_1dt \frac{1}{t^2 \sqrt{t^2-1}} e^{-\frac{2tr}{{\mkern0.75mu\mathchar '26\mkern -9.75mu\lambda}_e}}-1\right) \right]
\end{equation}
and contains the gradient of the scalar potential $\phi_{n;A}$ associated with nucleus $A$. The contribution from the electric form factor has been split into a high-frequency (HF) part
\begin{eqnarray}
\label{eqn:HF}
V^\textrm{HF}_A(r)= -A\left(Z_A,r\right)\frac{\alpha}{\pi} \phi_{n;A}(r) \int^{\infty}_1dt \frac{1}{\sqrt{t^2-1}}\left[ \left(1-\frac{1}{2t^2}\right)\right. \nonumber \\
\left.\times\left[\textrm{ln}\left(t^2-1\right)+4\textrm{ln}\left(\frac{1}{Z_A\alpha}+0.5\right)\right]-\frac{3}{2}+\frac{1}{t^2}\right]
e^{-\frac{2tr}{{\mkern0.75mu\mathchar '26\mkern -9.75mu\lambda}_e}},
\end{eqnarray}
obtained from theory after introducing a fixed low-frequency cutoff, and a low-frequency (LF) part
\begin{equation}\label{eqn:LF}
V^\textrm{LF}_A(r) = -B\left(Z_A\right) Z^4_A \alpha^5 m_e c^2 e^{-Z_Ar}.
\end{equation}
for which an approximate form has been chosen. Parameters $A\left(Z,r\right)$ and $B\left(Z\right)$ have been fitted to SE shifts of hydrogen-like ions \cite{Mohr_PhysRevA.45.2727,Mohr_PhysRevA.46.4421}.

These effective QED-potentials have been implemented by grafting routines from the numerical atomic code GRASP\cite{dyall1989grasp} onto the DFT-grid of the DIRAC code\cite{saue2020dirac,Saue_JCC2002,Sunaga2021qed}. Since the QED effects are principally generated in the vicinity of nuclei, on the order of reduced Compton wavelengths\cite{Artemyev_chap2016}, the numerical integration is only carried out on the nucleus associated with the QED-potential and with a radial cutoff. A problematic term in this respect is the low-frequency part of the Flambaum--Ginges potential since it by construction has the range of a hydrogen-like $1s$ orbital and for light elements may therefore overlap with other atoms in a molecule. Therefore, by default in the current implementation, QED-potentials are only activated for centers with $Z>18$.

The potentials can also be used at the eXact 2-Component relativistic (X2C) level\cite{ilias:jcp2007}. In the present work we employ the eXact 2-Component (X2C) molecular-mean field (mmf) Hamiltonian\cite{Sikkema_Visscher_Saue_Ilias_2009}. In this approach an initial HF calculation is carried out at the 4-component relativistic level. The converged Fock matrix is then transformed to 2-component form and employed in the ensuing correlated calculations.

\subsection{Parity-violating energy shift}\label{sec:PV}
Theoretical studies of parity violation in atoms and molecules are based on an effective Hamiltonian describing the weak interaction between electrons and nucleons mediated by the exchange of virtual $Z^0$ bosons. A detailed derivation and discussion of this Hamiltonian, starting at the level of quarks, is for instance found in Ref.~\cite{Greiner_weak} (the reader may also consult Refs.~\citenum{Berger:Parity-violation_effects_in_molecules,bast2011analysis,Bakasov_JCP_1998,halzen:quark,Hung_Sakurai_ARNPS1981,Commins_Bucksbaum_ARNPS1980}). The $Z^0$ boson is neutral, like its electroweak partner, the photon. On the other hand, it has mass --- 91.1876(21) GeV/c \cite{nakamura2010particle}, which is close to that of $^{98}\textrm{Mo}$, the most abundant isotope of molybdenum \cite{bast2011analysis}. This confers a short-range nature to the interaction. The effective Hamiltonian accordingly has the form of a contact interaction between two fermion 4-currents,
\begin{equation}
  \hat{h}_{\rm{eff}}\sim\int j^{\mu;e}\left(\boldsymbol{r},t\right)\left(\sum_{i}^{Z}j_{\mu;i}^{p}\left(\boldsymbol{r},t\right)+\sum_{i}^{N}j_{\mu;i}^{n}\left(\boldsymbol{r},t\right)\right)d^{3}\boldsymbol{r}
\end{equation}
in line with the form suggested by Fermi in 1934\cite{Fermi_ZfP1934}. The fermion currents of electromagnetic interaction have the form $\bar{\psi}\gamma^{\mu}\psi=\psi^{\dagger}(I_4,\boldsymbol{\alpha})\psi$, which is termed a vector current since the space-like part of the current transforms as a vector under space inversion. The fermion currents associated with the weak interaction are of mixed vector and axial vector (V-A) nature\cite{Sudarshan_PR1958,Feynman_PR1958,Sakurai_NC1958}, the latter on the form $\psi^{\dagger}(\gamma^5,\boldsymbol{\Sigma})\psi$, where appears
\begin{equation}
\gamma^{5}=\left[\begin{array}{cc}
0_{2} & I_{2}\\
I_{2} & 0_{2}
\end{array}\right].
\label{eqn:gamma5}
\end{equation}
For the nucleon currents a non-relativistic approximation is used, setting the small components of the 4-spinors to zero, which in practice means that terms involving the odd operators $\gamma^5$ and $\boldsymbol{\alpha}$ vanish whereas even operators are reduced to their 2-component counterparts, that is,  $I_4,\boldsymbol{\Sigma}\rightarrow I_2,\boldsymbol{\sigma}$.

Parity-violation arises from the $V_e-A_{p,n}$ and $A_e-V_{p,n}$ cross-terms. The former coupling leads to a nuclear spin-dependent effective Hamiltonian, which has so far been ignored in the calculation of parity-violation energies, but has been employed in studies of parity-violation in the NMR spectra of chiral molecules\cite{Barra_PLA1986,Barra_EPL1988,Laubender_CPC2003,Soncini_PhysRevA.68.033402,Bast2006JCP,Nahrwold_JCP2009}. Here we limit attention to the latter coupling, which leads to a nuclear spin-independent effective Hamiltonian on the form
\begin{equation}
\hat{h}^{\rm{PV}}_{\rm{eff}} = \frac{{{G_{\rm{F}}}}}{{{\rm{2}}\sqrt 2 }}\sum\limits_A {Q_{W,A}}\sum\limits_i {\gamma_i^5\rho_A\left({\mathbf{r}_i}\right)} .
\label{eqn:pv_op}\end{equation}
where
\begin{equation}
Q_{W,A} =  - N_A + Z_A\left( {1 - 4{{\sin }^2}{\theta_W }} \right),
\end{equation}
is the weak charge\cite{Bouchiat_JPF1974} of nucleus $A$ given in terms of the numbers $N_A$ and $Z_A$ of neutrons and protons, respectively, as well as the weak mixing angle $\theta_\textrm{W}$ \cite{glashow1961partial,weinberg1967model}. The current recommended value\footnote{The value depends on the renormalization scheme employed. The PDG employs the modified minimal subtraction scheme, whereas the conceptually simpler, but as such less accurate\cite{CODATA2002} CODATA value $\sin^2\theta_\textrm{W}=1-(m_W/m_Z)^2=0.22290(30)$ \cite{CODATA2018}, given explicitly in terms of the masses of the $Z^0$ and $W^{\pm}$ bosons is obtained from the on-shell scheme.} from the Particle Data Group (PDG)\cite{PDG2020} is $\sin^2\theta_\textrm{W}=0.231 22(4)$.

The presence of the Fermi coupling constant $G_\textrm{F}=2.22255\times 10^{-14}E_ha_0^3$ in the effective Hamiltonian \eqref{eqn:pv_op} shows that the interaction is truly weak. The operator can therefore not be added to the electronic Hamiltonian in an energy calculation, since its contribution would be masked by numerical noise, rather it will be obtained perturbatively as an expectation value. The expectation value with respect to the electronic wavefunction $\Psi$ can in turn be expressed as a sum over atomic contributions,
\begin{equation}
E^{\rm{PV}} = \sum\limits_A  E^{\rm{PV}}_A =  \frac{G_F}{2\sqrt 2}\sum\limits_A Q^A_{W} M_A^{\rm{PV}};
\quad
M_A^{\rm{PV}} = \sum\limits_i\left<\Psi \left| \gamma_i^5\rho_A\left({\mathbf{r}_i}\right) \right|\Psi\right>.
\end{equation}
The reduced PV matrix element $M_A^{\rm{PV}}$ is of highly local nature, not only because of the presence of  the normalized nuclear charge density $\rho_A$ of each center $A$, but also because the $\gamma_5$ matrix \eqref{eqn:gamma5} couples large and small components, the latter being of highly local atomic nature for positive-energy orbitals.

The parity-violating energy $E^{\rm{PV}}$ is strictly zero in the absence of spin-orbit coupling and its meaningful calculation in a non-relativistic framework hence requires double perturbation theory\cite{Gajzago_AIPConf1974,Rein_PLA1979,hegstrom1980calculation,mason1984parity,lazzeretti1997calculation,Bakasov_JCP_1998,berger2000multiconfiguration}. The non-relativistic form of the effective Hamiltonian \eqref{eqn:pv_op} and spin-orbit interaction then contributes factors $Z_A^3$ and $Z_B^2$, respectively, giving an overall $Z^5$  scaling. Single-center contributions ($A=B$) vanish in a non-relativistic framework\cite {hegstrom1980calculation}, but not necessarily in a relativistic setting\cite{bast2011analysis}.

In a previous study, we analyzed the parity-violating energy of the \ce{H2X2}-series of molecules at the Hartree--Fock(HF) and Kohn--Sham(KS) levels\cite{bast2011analysis}. The $E^{\rm{PV}}$ expectation value was decomposed in terms of pre-calculated atomic orbitals. It was
found that a contribution $E^{\rm{PV}}_A$ is completely dominated by contributions from the same center $A$. Inserting the form of relativistic
atomic orbitals
\begin{equation}
  \psi(\mathbf{r})=\left[
  \begin{array}{r}R^L(r)\chi_{\kappa,m_j}(\theta,\phi)\\iR^S(r)\chi_{-\kappa,m_j}(\theta,\phi)\end{array}
  \right],
\end{equation}
the expectation value $E^{\rm{PV}}_A$ can be expressed in terms of products of radial and angular contributions.
From the angular part of  it was found that only pairs of atomic orbitals with \textit{opposite} values of $\kappa$ and the \textit{same} value of $m_j$ can contribute. From the small $r$ behaviour of the radial functions it was found that for point charge nuclei only $(s_{1/2},p_{1/2})$ - pairs of atomic orbitals can contribute, whereas for extended nuclei contributions with $|\kappa|>1$ are allowed. Finally, for the \ce{H2X2}-series of molecules we found that the reduced PV matrix element $M^X_{\rm{PV}}$ can be remarkably well approximated by the mixing of valence $ns_{1/2}$ and $np_{1/2}$ orbitals on the same center $X$, that is
\begin{equation}
  M^X_{\rm{PV}}\approx M^X_{\rm{PV};\it{nn}}\times 2\rm{Re}\left[\sum_{\it{i}} \it{c}(\it{ns}_{\rm{1/2}}^X)^\ast_i \it{c}(\it{np}_{\rm{1/2}}^X)_i\right],\label{eq:MX_PV}
\end{equation}
where
\begin{equation}
M^X_{\rm{PV};\it{mn}}=\langle ms_{1/2}^X|\gamma_5\rho^X|np_{1/2}^X\rangle =  \langle R_{ms_{1/2}}^{L;X}|\rho^X|R_{np_{1/2}}^{S;X}\rangle_r + \langle R_{ms_{1/2}}^{S;X}|\rho^X|R_{np_{1/2}}^{L;X}\rangle_r .\label{eq:PVint}
\end{equation}
In \eqref{eq:MX_PV} $c(\psi)_i$ is the coefficient of atomic orbital $\psi$ in the expansion of molecular orbital $i$. The quantity $2\rm{Re}\left[\sum_{\it{i}} \it{c}(\it{ns}_{\rm{1/2}}^X)^\ast_i \it{c}(\it{np}_{\rm{1/2}}^X)_i\right]$ thereby determines the mixing of atomic orbitals $ns_{1/2}$ and $np_{1/2}$ \textit{contributing} to the PV energy, but it should be noted that this particular form refers to a specific phase convention for atomic orbitals. Further details are found in Ref.~\citenum{bast2011analysis}.

\section{Computational details}\label{sec:comp}
Unless otherwise stated all calculations have been carried out using the DIRAC program \cite{saue2020dirac} (revisions 9e12f6d89, 13d3f86a3, 9533dd59e, 09012db84, and ccb40d212).
We employed a Gaussian model for the nuclear charge distribution in all calculations\cite{visscher1997dirac}.
Molecular geometries were taken from Table I of Ref.~\citenum{laerdahl1999fully} and were optimized at the MP2 level using quasirelativistic pseudopotentials for the heavy elements.
In all calculations the dihedral angle was fixed at 45$^{\circ}$ of the P enantiomer. For the weak mixing angle, we have used the value $\sin^2\theta_\textrm{W}=0.2319$ in line with previous works\cite{Bakasov_JCP_1998,laerdahl1999fully} and corresponding to the PDG 1994 recommended value\cite{PDG_1994}.

Various flavours of Dyall basis sets\cite{dyall2002relativistic, dyall2006relativistic, dyall2016relativistic} have been employed for all elements of the target molecules. These basis sets have been extended by tight exponents for better description of the PV operator.
These optimizations were carried out using the simplex method\cite{pressnumerical}, minimizing the target function
\begin{equation}\label{eq:PVerr}
\Delta M_{\rm{PV}}\left( \%  \right) = 100 \times \frac{1}{N}\sum\limits_{nm} {{\rm{abs}}\left( \frac{M_{\rm{PV,}\it{nm} }^{{\rm{num}}} -  M_{\rm{PV,}\it{nm} }^{\rm{bas}}}{M_{\rm{PV,}\it{nm} }^{\rm{num}}} \right)} ,
\end{equation}
where $M_{\rm{PV,}nm }^{{\rm{num}}}$ and $M_{\rm{PV,}nm }^{\rm{bas}}$ refer to reduced PV matrix elements, \eqref{eq:PVint}, calculated numerically or in the current basis set, respectively.
For this, we employed a development version of the atomic GRASP code\cite{dyall1989grasp} allowing the use of Gaussian basis sets, in this case the Dyall.$n$zp ($n=2,3,4$) basis sets for SCF calculations.

For the MP2 and CCSD calculations a computational protocol was developed based on a calibration study of the parity-violating energy $E_{PV}$ of the \ce{H2Se2} molecule,
described in the Supplemental Material. All calculations employed the molecular-mean-field X2C Hamiltonian ($^2DC_M$) \cite{Sikkema_Visscher_Saue_Ilias_2009} based on the  Dirac--Coulomb  Hamiltonian. For technical reasons, MP2 and CCSD analytical expectation values were calculated using the RELCCSD\cite{Visscher_Lee_Dyall_1996,shee2016analytic} and ExaCorr\cite{pototschnig2021implementation} modules, respectively. The expectation
values were calculated using the Lagrangian-based analytical energy derivative technique described in Ref.~\citenum{shee2016analytic} and references therein.
For the CCSD calculations, we used the dyall.aae3z basis sets \cite{dyall2002relativistic, dyall2006relativistic, dyall2016relativistic}, extended with tight functions, as described above. We correlated $\textrm{(}n-1\textrm{)}s\textrm{(}n-1\textrm{)}p\textrm{(}n-1\textrm{)}dnsnp$ electrons, and truncated the virtual space at 100 $E_h$.
For MP2 the calibration study indicated that only valence orbitals need to be correlated, so we used the the dyall.acv3z basis.

For comparison we carried out DFT calculations using the extended dyall.aae3z basis sets and functionals from different rungs of the Jacob's ladder of density functional approximations\cite{perdew:jacob}: LDA (SVWN5)\cite{Slater1951LDA, Vosko1980LDA}, BLYP\cite{Becke1988BLYP, Lee1988BLYP, Miehlich1989BLYP}, and B3LYP\cite{Stephens_Devlin_Chabalowski_Frisch_1994,Becke1993B3LYP,Hertwig1997B3LYP}, which is Becke's series, and PBE\cite{Perdew1986PBE, Perdew1996PBE} and PBE0\cite{Adamo1999PBE0, Ernzerhof1999PBE0}, which is Perdew's series. We also used the long-range correlated functional CAMB3LYP\cite{Yanai2004}.

\section{Results and discussion.}\label{sec:results}
For convenience, we have collected our HF, MP2 and CC results in Table 6 of the Supplemental Material together with previous results\cite{Thyssen_PRL2000,van2005first,shee2016analytic} calculated at the same geometries\cite{laerdahl1999fully}. Our outputs are available via the zenodo repository\cite{paper:dataset}.

\subsection{Optimization of tight exponents}\label{Sec:simplex}
In Table \ref{Tab:error_M_PV} we show the average relative error $\Delta M_{\rm{PV}}$  \eqref{eq:PVerr} of reduced PV matrix elements for various basis sets with respect to the reference numerical calculation. For Se, we see a reduction from 19.8\% to 3.2\% when going from dyall.2zp to dyall.4zp unmodified basis sets. At the dyall.3zp level the relative error is reduced for the heavier elements since these basis sets already contain quite tight exponents.
Adding tight exponents, however, the error can be brought below 1\% at all basis set levels.
Specifically, we added $2p$ for O, $1s3p$ for S and $1p$ for Te at the dyall.3zp level. For the Se atom, we added $1s3p$ for dyall.2zp, $2p$ for dyall.3zp, and $1p$ for dyall.4zp, respectively. The values of the optimized exponents are given in the Supplemental Material. For Po no exponents were added since the relative error was already well below 1\%.

\begin{table}[]
\tbl{Average relative error  $\Delta M_{\rm{PV}}$(\%) in reduced PV matrix elements in different basis sets (with and without added tight exponents) with respect to the reference numerical values.}
{\begin{tabular}{lcccccccc}
\hline
 & O & S & \multicolumn{3}{c}{Se} & Te & Po\tabularnewline
basis & dyall.3zp & dyall.3zp & dyall.2zp & dyall.3zp & dyall.4zp & dyall.3zp & dyall.3zp\tabularnewline \hline
without & 11.4 & 11.4 & 19.8 & 8.5 & 3.2 & 3.4 & 0.1\tabularnewline
with & 0.9 & 0.5 & 0.9 & 0.6 & 0.2 & 0.1 & \tabularnewline
\hline
\end{tabular}}
\label{Tab:error_M_PV}
\end{table}

\subsection{Comparison of methods}
Parity-violating energies of the \ce{H2X2} molecules, calculated at the 45$^{\circ}$ dihedral angle of the P enantiomer using different methods, are  collected in Table \ref{Tab:DFT}. Using the CCSD results as reference, we find that the mean relative error (MRE) of the HF method, corresponding to the relative correlation effect, is -2.8\%, so rather modest, in line with previous studies\cite{Thyssen_PRL2000,shee2016analytic}. On the other hand, the standard deviation $\Delta_{\rm{std}}$ approaches ten percent, showing quite a spread. MP2 has a slightly smaller MRE of opposite sign (2.0 \%) and similar spread. No density functional outperforms HF; their average MRE and $\Delta_{\rm{std}}$ are -11.1 \% and 7.0\%, respectively. The best performing functional is PBE0 (-5.7$\pm$ 5.0 \%), the worst CAMB3LYP (-17.2$\pm$4.3 \%). These observations are as such disappointing since, as already mentioned, the computational protocol developed for the scan of candidate molecules for the Paris experiment is based on DFT. On the other hand, the protocol calculates a vibrational shift, corresponding to the PV energy difference of two vibrational levels, rather than an absolute electronic PV energy as here. Further study is therefore required to judge the performance of the computational protocol.

\begin{table}[h]
\tbl{Summary of the values of $E_{\rm{PV}} (E_h)$, calculated at the 45$^{\circ}$ dihedral angle of the P enatniomer of \ce{H2X2} molecules at various levels of theory. "err" is the percentage error with respect to CCSD. }
{\begin{tabular}{lrrrrrrrrrr}
\hline
 & \multicolumn{2}{c}{\ce{H2O2}} & \multicolumn{2}{c}{\ce{H2S2}} & \multicolumn{2}{c}{\ce{H2Se2}} & \multicolumn{2}{c}{\ce{H2Te2}} & \multicolumn{2}{c}{\ce{H2Po2}}\tabularnewline
 Method & $E_\textrm{PV}$ &  err (\%) & $E_\textrm{PV}$ & err (\%) & $E_\textrm{PV}$ & err (\%) & $E_\textrm{PV}$ & err (\%) & $E_\textrm{PV}$ & err (\%) \tabularnewline \hline
CCSD     & -5.917E-19 &       & -2.166E-17 &         & -2.463E-15 &       & -3.536E-14 &       & -1.337E-12 &      \tabularnewline
MP2      &  -6.525E-19 &  -10.3 & -2.067E-17 &  4.6 & -2.191E-15 & 11.1  & -3.236E-14 &   8.5 &  -1.390E-12&  -3.9\tabularnewline
HF       & -6.808E-19 & -15.0 &  -2.160E-17 & 0.3 & -2.301E-15 &  6.6  & -3.403E-14 &   3.8 & -1.467E-12 &   -9.7\tabularnewline
LDA      & -6.581E-19 & -11.2 & -2.659E-17 &  -22.7& -2.910E-15 & -18.2 & -4.164E-14 & -17.7 & -1.308E-12 &    2.2\tabularnewline
BLYP     & -6.303E-19 & -6.5  & -2.663E-17 & -22.9& -2.864E-15 & -16.3 & -4.095E-14 & -15.8 & -1.324E-12 &    1.0\tabularnewline
B3LYP    & -6.515E-19 & -10.1 & -2.576E-17 &  -18.9 & -2.796E-15 & -13.5 & -4.032E-14 & -14.0 & -1.372E-12 &   -2.6\tabularnewline
CAMB3LYP & -6.627E-19 & -12.0 & -2.621E-17 & -21.0& -2.941E-15 & -19.4 & -4.263E-14 & -20.5 & -1.515E-12 &  -13.3\tabularnewline
PBE      & -6.269E-19 & -5.9  & -2.498E-17 & -15.3 & -2.653E-15 & -7.7  & -3.814E-14 &  -7.8 & -1.261E-12 &    5.7\tabularnewline
PBE0     & -6.487E-19 & -9.6  & -2.404E-17 & -11.0 & -2.570E-15 & -4.4  & -3.729E-14 &  -5.4 & -1.313E-12 &    1.8\tabularnewline
\hline
\end{tabular}}
\label{Tab:DFT}
\end{table}

\subsection{Radiative corrections to PV energies}

In Table \ref{Tab:QED} we show QED-effects on parity-violating energies for the heavier members of the \ce{H2X2}-series, calculated at the CCSD level. We used the Uehling term (U) for the inclusion of vacuum polarization (VP) and the Pyykk{\"o}--Zhao (PZ) or Flambaum--Ginges (FG) potentials for electron self-energy (SE). Fitting our data to a simple power law $\sim Z_X^\gamma$ with respect to the nuclear charge of the heavy atom $X$, we find that the uncorrected PV energy has exponent 7.0, whereas the radiative corrections show stronger scaling, with $\gamma$ equal to 8.1 and 7.7, for $\Delta$(U+PZ), and $\Delta$(U+FG), respectively. In passing we note that this trend may not hold for superheavy systems, since SE and VP terms tend to have opposite sign, as known for valence atomic properties (see for instance Refs.~\citenum{pyykko1998estimated,Thierfelder2010PRA}).

From Table \ref{Tab:QED} we see that the radiative correction may reach a few percent of the PV energy, but only if we look at the results obtained with the FG potential. The discrepancy between the results obtained with the PZ and FG potentials could signal that one or both are being used beyond their domain of validity, since they have been principally designed to reproduce QED energy shifts.

\begin{table}[]
  \tbl{QED effects on the parity-violating energy $E_\textrm{PV}$ ($E_h$), calculated at the 45$^{\circ}$ dihedral angle of the P enantiomer of $\textrm{H}_2\textrm{X}_2$ molecules.
The radiative corrections have been included using the Uehling potential (U) for vacuum polarization and either the Pyykk\"{o}--Zhao (PZ) or Flambaum--Ginges (FG) effective potentials for electron self-energy. Absolute and relative differences are also shown.}{\begin{tabular}{lccccccc}
\hline
 & non-QED & U+PZ & U+FG & $\Delta$(U+PZ)($E_h$) & $\Delta$(U+FG)($E_h$) & $\Delta$(U+PZ)(\%) & $\Delta$(U+FG)(\%) \tabularnewline  \hline
\ce{H2Se2} & -2.463E-15 & -2.458E-15 & -2.434E-15 & 4.6E-18 & 2.9E-17 & -0.19 & -1.18\tabularnewline
\ce{H2Te2} & -3.536E-14 & -3.525E-14 & -3.471E-14 & 1.2E-16 & 6.6E-16 & -0.33 & -1.89\tabularnewline
\ce{H2Po2} & -1.337E-12 & -1.331E-12 & -1.306E-12 & 6.8E-15 & 3.1E-14 & -0.51 & -2.38\tabularnewline
\hline
\end{tabular}}
\label{Tab:QED}
\end{table}

\section{Conclusions}\label{sec:conc}
In this article we have given improved values for the parity violating energy of \ce{H2X2} molecules (X = O, S, Se, Te, Po). Improvements include: i) Property-optimized basis sets, which go beyond conventional energy-optimized basis sets. ii) We have used a new relativistic coupled cluster code, ExaCorr, geared towards massively parallel calculations\cite{pototschnig2021implementation}, allowing extensive calibration studies in order to elaborate a CCSD computational protocol for PV energies and to include \ce{H2Po2} among the target molecules iii) Inclusion of QED-effects through effective potentials.

Using our CCSD numbers as reference, we find the performance of DFT somewhat disappointing. This is possibly worrisome, since DFT is the basis of the computational protocol elaborated for the scan of candidate molecules for the high-resolution laser spectroscopy experiment which is under development in Paris and which aims to provide the first observation of parity violation in molecules as a shift in the vibrational spectra of enantiomers of chiral molecules. We shall pursue our investigation of the reliability of this protocol by calculating at the coupled cluster level the vibrational shift, involving energy differences, of candidate molecules. On the positive side, it has been noted that HF, here better performing, tends to predict larger PV vibrational shifts than DFT \cite{DeMontigny_PCCP2010}.

A major novelty of this work is the first calculation of molecular PV energies $E_{\rm{PV}}$ including QED effects. For the heaviest molecules in the \ce{H2X2}-series the radiative corrections are on the order of a few percent. However, we do note that this depends on what effective self-energy potential is used. This may indicate limitations to their applicability, an issue that encourages us to develop more rigid approaches to incorporate QED effects.

Further improvements would include going beyond the zeroth-order term, that is, the instantaneous Coulomb term, of the fully relativistic electron-electron interaction\cite{saue2011relativistic}.
The contribution from the Breit term, or more precisely spin-other-orbit interaction, has been estimated by Berger within the Breit--Pauli framework\cite{Berger_JCP2008}. The Breit effects for \ce{H2X2} (X = O, S, Se, Te, Po) are 12\%, 4.7\%, 1.9\%, 1.2\%, 0.7\%, respectively and hence monotonically decreasing along the series, contrary to the radiative corrections.

Another improvement would be to increase the excitation level beyond that of CCSD. An indication of the role of higher-order correlation effects comes from the finite-field study of Thyssen \textit{et al}\cite{Thyssen_PRL2000}: The relative difference between CCSD(T) and CCSD was found to be $-1.0 \%$ and $-2.5 \%$ for \ce{H2O2} and \ce{H2S2}, respectively. This suggests that these effects may be at the percentage level and increase down the series. This clearly calls for further study.

\section*{Acknowledgement(s)}
This paper is dedicated to John F. Stanton, a driving force of coupled cluster theory for 35 years !
We would like to thank David Sanchez(Toulouse) and Johann Pototschnig(Vienna) for a lot of technical support.
This work is part of the project molQED funded by the Agence Nationale de la Recherche (ANR, France).
This research used resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC05-00OR22725. We are thankful to CALMIP (Calcul en Midi-Pyr\'{e}n\'{e}es) and LCPQ (Laboratoire de Chimie et Physique Quantiques) clusters for the additional computation time. We are also thankful to the supercomputer of ACCMS (Kyoto University) and Research Institute for Information Technology, Kyushu University (General Projects) for test calculations.
A.S. acknowledges financial support from the Japan Society for the Promotion of Science (JSPS) KAKENHI Grant No. 17J02767 and 20K22553, and JSPS Overseas Challenge Program for Young Researchers, Grant No. 201880193.


\section{Introduction}\label{sec1}\blfootnote{This is a preprint manuscript accepted for publication at ECML-PKDD 2021.}

\noindent
The majority of financial exchange venues utilise a continuous double auction (CDA) mechanism \cite{friedman1993double} for matching orders. Under CDA formation, both {\em ask} orders (orders to sell a given quantity at a given price) and {\em bid} orders (orders to buy a given quantity at a given price) arrive at the venue continuously, with no minimum time interval limit. When a new order arrives, if it does not immediately execute, it will enter the limit order book (LOB); which contains a list of current bids and a list of current asks, both sorted by price-time priority. Therefore, the LOB contains valuable information on the instantaneous demand and supply for a particular financial asset (e.g., a stock, a commodity, a derivative, etc.). For this reason, LOB data has been used for many and various studies, including exploration of the price formation mechanism \cite{parlour1998price}, market anomaly detection \cite{ye2019extracting}, and testing of  trading algorithms \cite{abergel2020algorithmic}.

However, there remain some obstacles for the wider application of LOB data. Firstly, LOB data subscription fees are usually high, sometimes amounting to tens of thousands of dollars per annum.\footnote{http://www.nasdaqtrader.com/Trader.aspx?id=DPUSdata} This might be a trivial sum for an institutional subscriber, however for individual investors and researchers this significant expense can hold them back. Further, LOB data is entirely unavailable in venues that deliberately do not make order information public, for instance some e-commercial markets and {\em dark pools} (e.g., see \cite{MPC-CarSmaTal-19,EPRINT:CarSmaTal20}). This challenge attracts researchers to consider the possibility of recreating the LOB from a more easily available source, such as {\em trades and quotes} (TAQ) data. TAQ data contains the {\em top price level} information of a LOB (the lowest-priced ask and highest-priced bid), together with a history of transactions. It is published to the public for free in most venues. Blanchet et al. \cite{blanchet2017unraveling} have previously demonstrated that it is possible to predict daily average order volumes resting at different price levels of the LOB, using only TAQ data for parameter estimation. More recently, from a deep learning perspective, the LOB recreation model (LOBRM) was proposed to formalize the task as a time series prediction problem, and an ensembled recurrent neural network (RNN) model was successfully used to predict order volumes in a high frequency manner for the first time \cite{shi2021lob}. Nevertheless, there exist two key restrictions in the LOBRM study: (1) The original LOBRM study was conducted in an interpolation style on only one day's length of LOB data, for two stocks. For the model to be applied in a real world application scenario, such as online prediction of market price movements, LOBRM performance requires evaluation on an extended multi-day dataset, with chronological training and testing such that there is no possibility of lookahead bias; (2) The ordinary differential equation (ODE) kernel used in the original LOBRM model has high computation complexity and is therefore inefficient for more realistic application scenarios when large amounts (weeks or months) of training data is used.

\paragraph{\bf Contributions:}
\begin{enumerate}
\item We advance the workflow and structure of the LOBRM model, such that: (i) a time-weighted z-score standardization for LOB features is used to enhance the model's generalization ability; and (ii) the original ODE kernel is substituted for an exponential decay kernel to enable faster inference of latent states, greater runtime efficiency, and a reduction in overfitting.

\item We use chronological training and testing to conduct experiments on an extended LOBSTER dataset that is an order of magnitude larger than the original dataset. We find that: (i) LOBRM with continuous decay kernel is superior in modelling the irregularly sampled LOB; and (ii) the module ensembling of LOBRM is effective.
\item We draw new empirical findings that further enrich the current literature: (i) the proposed sparse encoding method for TAQ data has good generalization ability and can facilitate manifold tasks including LOB prediction and price trend prediction; (ii) prediction accuracy of the LOBRM is negatively related to volume volatility at unseen price levels; and (iii) the influence of stochastic drift on model performance can be alleviated by increasing the amount of historical training samples.
\end{enumerate}

\section{Background and Related Work}\label{sec2}
\noindent
\subsection{The Limit Order Book (LOB)}
In a CDA market, bids and asks with specified price and quantity (or {\em volume}) are submitted, cancelled, and transacted continuously. The
LOB contains an ask side and a bid side, with ask orders arranged in price ascending order and bid orders arranged in price descending order. Ask orders with the lowest price ({\em best ask}) and bid orders with the highest price ({\em best bid}) form the top level of a LOB, and their respective prices are called {\em quotes}. If a newly submitted ask (or bid) price is not higher (or lower) than the best bid (or ask), a trade happens. TAQ data contains all historical quotes and trades in the venue. That is, LOB data contains strictly more information than TAQ data. Fig.~\ref{fig:LOB} provides a visual illustration of a LOB and the relationship between the LOB and the TAQ data.

\begin{figure}[t]
  \centering
  \includegraphics[width=0.85\linewidth]{fig_lob.pdf}
  \caption{A LOB of four price levels evolving with time. White and blue boxes indicate the top level and deeper levels of the LOB. Grey boxes indicate market events stream. Orange box indicates trade records. White and orange boxes together form TAQ.}
  \label{fig:LOB}
\end{figure}

Traditional statistical models of the LOB assume that LOB evolution follows the rules of a Markovian system, with market events (order submission, cancellation, and transaction) following stochastic point processes, such as a Hawkes process or a Poisson process \cite{abergel2015long,cont2013price}. This formulation generalizes a LOB market of high complexity to a dynamic system controlled by a few parameters, where probabilistic theorems like the law of large numbers \cite{horst2017weak} and stationary equilibrium in a Markovian system \cite{blanchet2017unraveling} can be utilized to draw long-term empirical conclusions. However, while statistical modelling can capture long term behaviour patterns of the LOB, these models cannot consistently perform well in a high frequency domain.
In recent years, there has been an emergence of research using deep learning approaches to model and exploit the LOB. Sirignano et al. \cite{sirignano2019universal} performed a significant study on a comprehensive pool of 500 stocks. They revealed that features learned by a Long Short-Term Memory (LSTM) network can be utilized to predict next mid-price movement direction $\{$\textit{up, down}$\}$ with accuracy range $[0.65,0.76]$ across all 500 stocks. Their study also demonstrated that deep learning models suffer less from problems such as stochastic drift that exist in statistical models of the LOB. Other deep learning studies of the LOB include extracting high frequency indicators \cite{passalis2020temporal,tsantekidis2020using}, predicting future stock price movement \cite{makinen2019forecasting,zhang2019deeplob}, and training reinforcement trading agents \cite{kumar2020deep,nevmyvaka2006reinforcement}.

\subsection{Generating Synthetic LOB Data}
\noindent
Synthetic LOB data, generated by models that learn from the real LOB or imitate the {\em stylized facts} of a CDA market, has been used as an alternative when real LOB data is unavailable. The advantages of using synthetic LOB data lie in its low cost and infinite availability. It has been widely adopted to backtest trading algorithms, explore market dynamics, and facilitate teaching activities.

Synthetic LOB data can be generated using three mainstream methodologies.
\textbf{Agent-based models} have been well studied and are the most popular approach for generating a synthetic LOB.
By configuring agents that trade using common strategies, such as market makers, momentum traders, and mean reversion traders, the synthetic LOB can closely approximate the stylised facts of a real LOB \cite{mcgroarty2019high}. Stock market simulators have a long history, from the Santa Fe artificial stock exchange \cite{arthur1996asset} to recent multi-agent exchange environments \cite{belcak2020fast}.
\textbf{Generative models} attempt to learn regularities embedded in market event streams or the LOB directly. One representative research by Li et al. \cite{li2020generating} utilized generative adversarial networks to learn and replicate the historical dependency among orders. The synthesised order stream and resulting LOB were found to  closely resemble the real market data.

We consider the aforementioned two approaches as {\em unsupervised}, since no real LOB data is used to verify the authenticity of the generated data. Model quality can only be verified by testing whether certain {\em stylized facts} exist in the synthetic data.
In contrast, \textbf{supervised models} use real LOB data as ground truth. As indicated in \cite{blanchet2017unraveling}, TAQ data is informative of LOB volumes for small-tick stocks. By modelling the tail probability of price change per trade in a Markovian LOB, daily average order volumes at unseen price levels can be estimated by the steady-state distribution of the infinite server queue. Shi et al. \cite{shi2021lob} further formulate the task of generating a synthetic LOB as a time series prediction problem using a continuous RNN. The proposed model (LOBRM) is able to predict LOB order volumes using a defined length of TAQ data as input. As long as a historical TAQ trajectory is available, the model is able to produce a historical replay of the LOB based on the knowledge it learned from supervised training. This paper concentrates on a further exploration of the LOBRM model presented in \cite{shi2021lob}.

\section{Model Formulation}\label{sec3}
\subsection{Motivation}\label{sec:motivation}
\noindent
The LOBRM model represents the first attempt to synthesize LOB data from a supervised deep learning perspective \cite{shi2021lob}. LOBRM is essentially an ensemble of RNNs that take TAQ data as input and produce LOB volume predictions as output. However, in the original study, there were three restrictions present:
    (1)
    Experiments were performed using a relatively small LOB dataset consisting of only one day's LOB data for two small-tick stocks. To verify its generalization ability, LOBRM requires testing on multi-day LOB data for a variety of stocks.
    (2) Experiments adopted a non-chronological approach to the formation of time series samples, such that samples were shuffled before splitting into training and testing sets. We intend to test model performance using a strictly chronological approach to ensure that LOBRM is applicable to real world online scenarios, with no possibility of introducing lookahead bias. Specifically, we use the first three days' data for training, the fourth day's data for validation, and the fifth day's data for testing. Time series samples are not shuffled, thus ensuring that chronological ordering is preserved.
    (3) The core module of the original LOBRM made use of an ODE-RNN, a RNN variant with ODE kernels to derive fine-grained time-continuous latent states \cite{rubanova2019latent}. As the computational complexity of ODE-RNN is $n$ times that of a vanilla RNN -- where $n$ is the granularity of the latent state -- it is not efficient for training data of large size. Therefore, we substitute the ODE-RNN for an RNN-decay module \cite{che2018recurrent}, which has been shown to be a more time-efficient model for irregularly sampled time series \cite{lechner2020learning}.

\subsection{Problem Description}
\noindent
To simplify the problem of recreating the LOB, we make the same assumptions proposed in \cite{shi2021lob}: (1) Following common practice (e.g., \cite{blanchet2017unraveling,horst2017weak}), the bid and ask sides of the LOB are modelled separately; (2) We only consider instantaneous LOB data at the time of each trade event, and ignore the LOB for all other events, such as order submission and cancellation; (3) We only consider the top five price levels of the LOB; (4) We assume that the price interval between each price level is exactly one {\em tick} -- the smallest increment permitted in quoting or trading a security at a particular exchange venue -- which is supported by empirical evidence that orders in the LOB for small-tick stocks tend to be densely distributed around the top price levels \cite{blanchet2017unraveling}. Following these assumptions, the price at different LOB levels can be directly deduced from the known quote price at target time. Therefore, the LOB recreation task resolves to the simpler problem of only predicting order volumes resting at each price level.

For generalization, we denote trades and quotes streams as $\left \{ TD_{i} \right \}_{i\in \textsl{n}}$ and $\left \{ QT_{i} \right \}_{i\in \textsl{n}}$ respectively,
and trajectories of time points for TAQ records as $\left \{ T_{i} \right \}_{i\in \textsl{n}}$ indexed by $\textsl{n}=\left \{ 1,\ldots,N \right \}$, where $N$ equals the number of time points in TAQ.
The LOB sampled at $\left \{ T_{i} \right \}_{i\in \textsl{n}}$ are denoted as $\left \{ LOB_{i} \right \}_{i\in \textsl{n}}$. For each record at time $T_{i}$, $QT_{i}=\left( p_{i}^{a(1)}, v_{i}^{a(1)}, p_{i}^{b(1)}, v_{i}^{b(1)}\right)$, where $p_{i}^{a(1)}, v_{i}^{a(1)}, p_{i}^{b(1)}, v_{i}^{b(1)}$ denote best ask price, order volume at best ask, best bid price, and order volume at best bid, respectively.
$TD_{i}=\left(p_{i}^{td}, v_{i}^{td}, d_{i}^{td}\right)$, where $p_{i}^{td}, v_{i}^{td}, d_{i}^{td}$ denote price, volume, and direction of the trade, with $+1$ and $-1$ indicating orders being sold or bought.
$LOB_{i}=(p_{i}^{a(l)}, v_{i}^{a(l)}, p_{i}^{b(l)}, v_{i}^{b(l)})$ depicts the price and volume information at all price levels, with $l\in(1,...,L)$, here $L=5$.
From the aforementioned model assumptions, we have $p_{i}^{a(l)}=p_{i}^{a(1)}+(l-1)\tau$ and $p_{i}^{b(l)}=p_{i}^{b(1)}-(l-1)\tau$, where $\tau$ is the minimum tick size (1 cent in the US market). For a single sample, the model predicts $(v_{I}^{a(2)},...,v_{I}^{a(L)})$ and $(v_{I}^{b(2)},...,v_{I}^{b(L)})$ conditioned on the observations of $\left \{ QT_{i} \right \}_{{I-S}:I}$ and $\left \{ TD_{i} \right \}_{{I-S}:I}$, with $S$ being the time series sample size, i.e., the maximum number of time steps that the model looks back in TAQ data history.

\subsection{Formalized Workflow of LOBRM} \label{subsec:model}
\noindent
\subsubsection{LOB Data Standardization}
As we intend to apply the LOBRM model on LOBs of five days' length for different financial assets, data standardization is necessary for the model's understanding of data of various numerical scales. We perform time-weighted z-score standardization on all LOB volumes, based on the fact that the LOB is a continuous dynamic system with uneven time intervals between updates.
We use $\left \{v_{n}^{(l)}\right \}_{0:N}$ and $\left \{T_{n}\right \}_{0:N}$ to indicate volume trajectory on price level $l$, and affiliated timestamps $N$ as the total number of LOB updates for training. Time-weighted mean and standard deviation are calculated as:
\begin{equation}
    \mathrm{Mean}^{(l)} = \begin{matrix}\sum_{i=0}^{N-1}{\left ( T_{i+1}-T_{i} \right ) v_{i}^{(l)}} \end{matrix}/\left (T_{N}-T_{0}\right )\\
\end{equation}
\begin{equation}
    \mathrm{Std}^{(l)} = \left ( \begin{matrix}\sum_{i=0}^{N-1}{\left ( {T_{i+1}-T_{i}} \right )} \left ( v_{i}^{(l)}-\mathrm{Mean}^{(l)} \right )^{2} \end{matrix}/\left (T_{N}-T_{0}\right )  \right )^{\frac{1}{2}} \\
\end{equation}

From empirical observation, we witness that the volume statistics on deeper price levels $\left \{v_{n}^{(2-5)}\right \}_{0:N}$ have similar patterns, while those statistics deviate from volume statistics on the top price level $\left \{v_{n}^{(1)}\right \}_{0:N}$. In particular, the top level tends to have much lower mean and much higher standard deviation (e.g., see later  Table~\ref{tab:dataset statistics}). Thus, we treat `top' and `deeper' levels as two separate sets to standardize. As trade volumes $\left \{v_{n}^{tq}\right \}_{0:N}$ are discrete events and do not persist in time, we use a normal z-score standardization for trade data. To avoid lookahead bias, statistics are calculated without considering test data. Finally, using these statistics, LOBs for training, validation, and testing are standardized.

\subsubsection{Sparse Encoding for TAQ}
TAQ data contains multi-modal information, including order type ({\em bid} or {\em ask}), price, and volume. While under the formulation of LOBRM, only order volumes at derived price levels (i.e., {\em deeper} levels 2-5) are predicted. We use a one-hot positional encoding, such that only volume information is encoded explicitly; while price is indicated by the position of non-zero elements in the one-hot vector.

Take the encoding of an ask quote as an example. Conditioned on current best ask price $p_{I}^{a(1)}$ and best bid price $p_{I}^{b(1)}$ at $T_{I}$, we represent the ask quote record $(p_{I-s}^{a(1)},v_{I-s}^{a(1)})$ and bid quote record $(p_{I-s}^{b(1)},v_{I-s}^{b(1)})$ at $T_{I-s}$, $s \in \left \{0,\ldots,S\right \}$ as:
\begin{equation}
    \left\{
    \begin{array}{cl}
        O_{2k-1}^{aq}, \mathrm{where}\,o_{k+sp_{s}^{a}}=v_{I-s}^{a(1)} \\
        O_{2k-1}^{bq}, \mathrm{where}\,o_{k+sp_{s}^{b}}=v_{I-s}^{b(1)}
    \end{array} \right.
    \label{eqt_5}
\end{equation}
where $k\in R$, $sp_{s}^{a}=(p_{I-s}^{a(1)}-p_{I}^{a(1)})/\tau$ and $sp_{s}^{b}=(p_{I-s}^{b(1)}-p_{I}^{b(1)})/\tau$. $O_{2k-1}$ is a one-hot vector with dimension $1\times (2k-1)$; and $o_{sp}$ denotes the ${sp}$-th element of the vector.
The value of $k$ is chosen to cover more than 90\% of past quote price fluctuations, relative to the current quote price. Here $k=8$, which means historical quotes with relative price $[-7,+7]$ ticks are encoded into feature vectors.
Then, a trade record $(p_{I-s}^{td}, v_{I-s}^{td},  d_{I-s}^{td})$, is represented as:
\begin{equation}
    \left\{
    \begin{array}{cl}
        O_{2k-1}^{td}, \mathrm{where}\,o_{p_{I-s}^{td}-p_{I}^{a(1)}}=v_{I-s}^{td} & \mathrm{if} \,d_{I-s}^{td}<0 \\
        O_{2k-1}^{td}, \mathrm{where}\,o_{p_{I-s}^{td}-p_{I}^{b(1)}}=v_{I-s}^{td} & \mathrm{if} \,d_{I-s}^{td}>0 \\
    \end{array} \right.
    \label{eqt_7}
\end{equation}

Finally, those three features are concatenated into $(O^{aq},O^{bq},O^{td})$ and are used as input. It can also be a concatenation of four features, with ask and bid trade represented separately. Later, in experiment section~\ref{sec:encoding}, we show this sparse encoding method can achieve enhanced robustness in LOB volume prediction and price trend prediction.

\subsubsection{Market Event Simulator Module (ES)}
The ES module models the overall net order arrivals as inhomogeneous poisson processes, and predicts LOB volumes from a dynamic perspective. If an RNN structure is used to iteratively receive encoded LOB features at every timestep, its latent state can be deemed as reflective of market microstructure condition over a short historical time window. A multi-layer perceptron (MLP) layer can then be used to decode latent states directly into vectors representing net order arrival rates at each price level.

The original LOBRM model uses ODE-RNN, a continuous RNN variant that learns fine-grained latent state between discrete inputs in a data-driven manner, to model the continuous evolution of market conditions. In ODE-RNN, a hidden state $h(t)$ is defined as a solution to an ODE initial value problem. The latent state between two inputs can then be derived using an ODE solver as:
\begin{equation}
    \begin{array}{cl}
    \frac{dh(t)}{dt}=f_{\theta }(h(t),t) \; \mathrm{where} \; h(t_{0})=h_{0}
    \end{array}
    \label{eqt_1}
\end{equation}
\begin{equation}
    \begin{array}{cl}
        {h_{i}}'=\mathrm{ODEsolver}(f_{\theta}, h_{i-1}, (t_{i-1}, t_{i}))
    \end{array}
    \label{eqt_2}
\end{equation}
in which function $f_{\theta }$ is a separate neural network parameterized by $\theta$.

Even though ODE-RNN contributes most to prediction accuracy, it is of high computation complexity (see Section~\ref{sec:motivation}).
Therefore, for an efficient use of LOBRM, especially when trained on large amounts of data, or used for online prediction, the ODE-RNN is unsuitable. Also, we find in chronological experiments that the ODE-RNN tends to cause overfitting due to the fully flexible latent states.

Faced with these challenges, we propose the use of a pre-defined exponential decay kernel \cite{che2018recurrent,lechner2020learning}, instead of the ODE kernel in the ES module. The inference of latent states between discrete inputs is denoted as:
\begin{equation}
    {h_{i}}'= h_{i-1}*exp\left ( -f_{\theta}\left ( h_{i-1} \right ) * \Phi \left ( t_{i-1},t_{i} \right ) \right )
\end{equation}
\begin{equation}
    h_{i} = \mathrm{GRUunit}\left ({h_{i}}',x_{i} \right )
\end{equation}
where $f_{\theta}$ is a separate neural network parameterized by $\theta$, and $\Phi$ is a smoothing function for time intervals to avoid gradient diminishing. A GRU unit \cite{cho2014learning} is then used for instant updating of the latent state at input timesteps. The advantages of employing an exponential decay kernel are threefold: (1) It allows for efficient inference of latent states; (2) It is a continuous RNN that imitates continuity of market evolution and includes temporal information within the model structure itself, sharing the advantages of continuous RNN in modelling irregularly sampled time series; (3) It is less likely to cause overfitting as the kernel form is predefined, whereas the latent state evolution in ODE-RNN is fully flexible.

We derive the vector of net order arrival rates $\Lambda_{I-s}=[\lambda_{I-s}^{a(2)},...,\lambda_{I-s}^{a(L)} ]$ at time $T_{I-s}$ directly from the latent state $h_{I-s}$, using an MLP layer as:
\begin{equation}
    \Lambda_{I-s}=\mathrm{MLP} \left (h_{I-s} \right )
\end{equation}
After acquiring the trajectory of $\Lambda$ at all trade times over the defined length of time steps, we calculate the accumulated order volumes between $[T_{I-S},T_{I}]$ as:
\begin{equation}
    \begin{matrix} \sum_{i=I-S}^{I}\Lambda_i\times \Phi \left (T_{i+1}-T_{i} \right ) \end{matrix}
\end{equation}

\subsubsection{History Compiler Module (HC)}
The HC module predicts LOB volumes from a historical perspective. It concentrates on historical quotes that are most relevant to current LOB volumes at deeper price levels. More precisely, for an ask side model the volumes to be predicted at target time are of prices $\{p_{I}^{a(1)}+\tau, \ldots, p_{I}^{a(1)}+(L-1)\tau\}$. Only historical ask quotes with price within this range are used as inputs into the HC module.
As one-hot encoded ask quotes fall within the price range of $\{p_{I}^{a(1)}-k\tau, \ldots, p_{I}^{a(1)}+k\tau\}$, vectors need to be trimmed to remove verbose information and retain the most relevant information. Formally, we represent a trimmed ask quote record $(p_{I-s}^{a(1)},v_{I-s}^{a(1)})$ as:
\begin{equation}
    \left\{
    \begin{array}{cl}
    O_{L-1}, \mathrm{where}\,o_{sp_{s}^{a}}=v_{I-s}^{a(1)} & \mathrm{if} \,sp_{s}^{a}\in[1,L-1] \\ Z_{L-1} & \mathrm{otherwise}
    \end{array} \right.
    \label{eqt_3}
\end{equation}
where $O_{L-1}$ is a one-hot vector of dimension $1\times (L-1)$; $o_{sp_{s}^{a}}$ denotes the
${sp_{s}^{a}}$-th element of the vector; and  $Z_{L-1}$ denotes a zero vector of the same dimension. The intention of the HC module is to look back in history to check how many orders were resting at price levels we are interested in at target time. Thus, we manually trim input feature vectors to leave out verbose information. A discrete GRU unit is used to compile trimmed features and generate volume predictions which are used as supplements to ES predictions.

\subsubsection{Weighting Scheme Module (WS)}
This module is designed to combine the predictions from the ES and HC modules into a final prediction.
We follow the intuition that if the HC prediction for a particular price level is reliable from a historical perspective, a higher weight will be allocated to it, and {\em vice versa}. If quote history for a target time price level is both abundant and recent, we weight the information  provided on current LOB volumes as more reliable. The abundancy and timing of historical quotes are denoted by the masking sequence of HC inputs. Formally, we represent the mask of an ask side HC input as:
\begin{equation}
    \left\{
    \begin{array}{cl}
    O_{L-1}, \mathrm{where}\,o_{sp_{s}^{a}}=1 & \mathrm{if} \,sp_{s}^{a}\in[1,L-1] \\ Z_{L-1} &  \mathrm{otherwise}
    \end{array} \right.
    \label{eqt_4}
\end{equation}
A GRU unit is used to receive the whole masking sequence to generate a weighting vector of size $1\times\left ( L-1\right )$ to combine predictions from ES and HC modules.

\section{Experiment and Empirical Analysis}\label{sec4}
In this section, experiments are conducted on the extended LOBSTER dataset. These data were kindly provided by \url{lobsterdata.com} for academic research. The stocks and time periods were not selected by the authors. The dataset contains LOB data of five days' length for three small-tick stocks (Microsoft, symbol {\tt MSFT}; Intel, {\tt INTC}; and JPMorgan, {\tt JPM}). To ensure that there is no selection bias or ``cherry-picking'' of data, all three available small-tick stocks were used in this study.  The extended dataset is  approximately ten times the size of the dataset used in \cite{shi2021lob} and is a strict superset, therefore enabling easier results comparison.

\paragraph{Model specification:} The {\bf ES module} consists of a GRU with 64 units, a two layered MLP with 32 units and ReLU activation for deriving parameters of the decay kernel, and a two layered MLP decoder with 64 units and Tanh activation. The {\bf HC module} consists of a GRU  with 64 units, and a two layered MLP decoder with 64 units and LeakyReLU activation. The latent state size is set at 32 in ES and HC modules. The {\bf WS module} consists of a GRU with 16 units, and a one-layered MLP decoder with 16 units and Sigmoid activation. The latent state size is set at 16 in the WS module. L1 loss is used as the loss function; and all models are trained for 150 iterations, with a learning rate of 2e-4. Model parameters are chosen based on the lowest loss on the validation set.

\subsection{Data Preprocessing}
We clean the dataset by retaining only LOB updates at trade times, and removing LOB data during the first half-hour after market open and the last half-hour before market close as these periods tend to be volatile.
To alleviate the effect of outliers, we divide all volume numbers by 100 and winsorize the data by the range $[0.005,0.995]$. Then, the standardization method proposed in Section~\ref{subsec:model} is applied to the cleansed LOB data. Data statistics are illustrated in Table~\ref{tab:dataset statistics}.

\begin{table}[t]
\footnotesize
\centering
\caption{Volume statistics before standardization, showing time-weighted mean volume and standard deviation on the top level and deep levels (levels 2 to 5).}\vskip -0.5em
\label{tab:dataset statistics}
\setlength\tabcolsep{1.5mm}
\begin{tabular}{l cc cc cc cc}
\toprule
        & \multicolumn{2}{c}{MSFT}
        & \multicolumn{2}{c}{INTC}
        & \multicolumn{2}{c}{JPM}\\
        \rule{0pt}{2ex}
        & \multicolumn{1}{c}{Bid} & \multicolumn{1}{c}{Ask}
        & \multicolumn{1}{c}{Bid} & \multicolumn{1}{c}{Ask}
        & \multicolumn{1}{c}{Bid} & \multicolumn{1}{c}{Ask}\\
\midrule
Top     & \multicolumn{2}{c}{119.7/89.1} & \multicolumn{2}{c}{127.1/98.3} &  \multicolumn{2}{c}{32.7/38.5}\\
Deep   &  189.6/51.7   &  192.7/53.4  &  185.6/62.3  &  178.9/51.4  &  61.3/23.9  &  63.4/40.5 \\
\bottomrule
\end{tabular}
\end{table}

We use the first three days' data for training, the fourth day's data for validation, and the fifth day's data for testing. This is a significantly different approach to that taken in \cite{shi2021lob}, in which the task was essentially interpolation, as all time series samples were shuffled before splitting into training and testing sets (i.e., training and testing sets were not ordered chronologically). We extract TAQ data and labels directly from the standardized LOB. These data are then converted to time series samples using a rolling window of size $S$, such that the first sample consists of TAQ histories at timesteps $1$ to $S$ and is labeled by LOB volume at deep price levels at time step $S$. The second sample consists of TAQ history at timesteps $2$ to $S+1$ and is labeled using LOB at time step $S+1$, {\em et cetera}.
We set parameters $S=100$ and $k=8$, as in the original LOBRM model \cite{shi2021lob}.

\subsection{Model Comparison}
\noindent
In this section, we illustrate training results on generating synthetic LOB using mainstream regression and machine learning methods: (1) Support Vector Machine Regression with linear kernel (LSVR); (2) Ridge Regression (RR); (3) Single Layer Feedforward Network (SLFN); and (4) XGBoost Regression (XGBR). We then evaluate the performance of LOBRM with either discrete RNN or Continuous RNN module: (1) GRU; (2) GRU-T, with time concatenated input; (3) Decay; and (4) Decay-T, with time concatenated input.
Two criteria are used: (1) L1 loss on test set. As all labels are standardized into z-score, the loss indicates the multiple of standard deviations between prediction and ground truth; (2) R-squared, calculated on the test data using the method presented by Blanchet et al. \cite{blanchet2017unraveling}, to enable us to perform a strict comparison with the existing literature.

\begin{table}[tb]
\small
\centering
\caption{Model comparison. Criteria shown in format: test loss/R-squared; all numbers are in 1e-1. The lowest test loss and highest R-squared in each set are underlined.}\vskip -0.5em
\setlength\tabcolsep{0.8mm}
\begin{tabular}{l cc cc cc}
\toprule
        Model & \multicolumn{2}{c}{MSFT}
        & \multicolumn{2}{c}{INTC}
        & \multicolumn{2}{c}{JPM}\\
        \rule{0pt}{2ex}
        & \multicolumn{1}{c}{Bid} & \multicolumn{1}{c}{Ask}
        & \multicolumn{1}{c}{Bid} & \multicolumn{1}{c}{Ask}
        & \multicolumn{1}{c}{Bid} & \multicolumn{1}{c}{Ask}\\
\midrule
LSVR     &  10.00/0.29  &  9.74/0.36    &  8.72/0.21  &  10.30/0.14 &   6.91/0.54   &   3.96/0.51 \\
RR       &  7.90/0.55   &  7.70/0.67    &  6.18/0.53  &  7.00/0.40  &   6.86/0.65   &   4.10/0.46 \\
SLFN     &  7.06/0.68   &  6.88/0.88    &  5.55/0.86  &  5.88/0.53  &   6.61/0.71   &   3.69/0.62 \\
XGBR     &  6.50/0.43   &  6.91/0.80    &  \underline{5.21}/0.82  &  5.93/0.83  &   6.81/0.76    &   3.73/0.69 \\
LOBRM (GRU)    &  6.44/1.24   &  6.44/1.55    &  5.34/1.05  &  5.70/0.53  &   6.27/1.17  &   3.30/1.25  \\
LOBRM (GRU-T)  &  6.45/1.38   &  6.45/1.56    &  5.36/1.06  &  5.65/\underline{0.84}  &   \underline{6.13}/1.48  &   3.34/\underline{1.73}  \\
LOBRM (Decay)  &  6.54/1.30   &  6.47/1.49    &  5.33/\underline{1.11}  &  \underline{5.52}/0.73  &   6.35/1.20  &   \underline{3.29}/1.37  \\
LOBRM (Decay-T)&  \underline{6.28}/\underline{1.58}   &  \underline{6.15}/\underline{1.77}    &  5.33/1.07  &  5.64/0.82  &   6.18/\underline{1.54}  &   3.32/1.65  \\
\bottomrule
\end{tabular}
\label{tab:continuous}
\end{table}

\begin{figure}[t]
\centering
\begin{minipage}[t]{0.45\textwidth}
\centering
\includegraphics[width=0.99\textwidth]{loss_bid_plot.pdf}
\vskip -0.5em
\caption{Training (dash) and validation (line) loss curves for bid-side models.}
\label{fig:loss}
\end{minipage}
\hspace{1em
\begin{minipage}[t]{0.45\textwidth}
\centering
\includegraphics[width=0.99\textwidth]{corr.pdf}
\vskip -0.5em
\caption{Hourly test loss plotted against hourly volume standard deviation.}
\label{fig:corr}
\end{minipage}
\end{figure}

Table~\ref{tab:continuous} presents evaluation results. Judging from both criteria, we observe that non-linear models outperform linear models (LSVR and RR). The LOBRM family also outperforms traditional non-linear models (especially in terms of R-squared) by effective ensembling of RNNs.
This indicates that the recurrent structure of RNNs can facilitate the model's capability in explaining temporal variance in LOB volume.
Further, LOBRMs with a  continuous RNN module exhibit superior performance over those with a discrete RNN module in 4 out of 6 experiments. LOBRM (Decay-T) achieves the lowest average L1 test loss and highest average R-squared, indicating that incorporating temporal information in both feature vectors and latent state dynamics is the most suitable approach to capture the model's dependence on time.
Meanwhile we find the model tends to be overfitting, with validation loss starting to rise before training loss converges. Fig.~\ref{fig:loss} shows loss curves for LOBRM (Decay-T) trained with a higher learning rate of 2e-3. The phenomenon of overfitting also explains why we use a pre-defined exponential decay kernel instead of a fully flexible ODE kernel involving a lot more parameters to model temporal dependence.
Thus, in the following experiments, we continue using LOBRM (Decay-T) as the main model.

As shown in Table~\ref{tab:continuous}, test losses across 6 sets of experiments fluctuate in the range $[0.332,0.628]$ for LOBRM (Decay-T). To test whether this fluctuation results from order volume volatility, we use the pre-trained model to calculate the correlation of (hourly) test L1 loss across all experiments against the (hourly) standard deviation of non-standardized volume resting at deep price levels (see Fig.~\ref{fig:corr}).
We find that the correlation between hourly L1 loss and volume volatility is significant (p-value\textless0.01) with a correlation coefficient of $\rho=0.48$. Thus, we conclude that LOB prediction accuracy is negatively related to volume volatility.

Compared with the original LOBRM study, which uses a non-chronological training and testing method, we find that the performance of LOBRM trained in a chronological manner is weakened. In order to make the criteria comparable, we first transform the L1 loss in z-score into L1 loss on volume as a percentage of average volume, using $(zscore*std_{i})/mean_{i}$ for each dataset $i$.
The transformed loss on LOBRM (Decay-T) ranges in $[16.1\%,24.1\%]$ with a mean value of 18.9\%, whereas this criteria in non-chronological experiments has an average of 6.9\%. The main reason lies in that, in non-chronological experiments the model attempts to predict the volume at a target timestep when both volumes at past and future timesteps are known and used for training.
This leakage of information from the future introduces a form of lookahead bias into the model that tends to increase the prediction accuracy. In the experiments presented here, no future labels are used for training,  thus lookahead bias is eliminated.

We further compare our results of R-squared with results from the statistical model approaching a similar task proposed in \cite{blanchet2017unraveling}, in which presented R-squared values ranged in $[0.81,0.88]$ for daily average volume at different price levels. However, in comparison, our R-squared values presented in Table~\ref{tab:continuous} are regressed against the ground truth with no averaging procedure. As we have only one-day length of data for training, we cannot generate daily average volume prediction for R-squared calculation. When we average volume predictions into hourly frequency, the value of R-squared across 6 experiments ranges in $[0.48,0.88]$ with an average value of 0.71. We find our results comparable to the existing literature, even though those results are not measured using identical conditions.

\subsection{Ablation Study}
\noindent
We conduct an ablation study to demonstrate the effectiveness of module ensembling. Four experiments are conducted: HC, using only the history compiler; ES, using only the market events simulator; HC+ES, which includes the history compiler and event simulator, using a pre-defined weight to combine outputs; and HC+ES+WS, which is the full LOBRM with adaptive weighting scheme.

\begin{table*}[t]
\small
\centering
\caption{Ablation study: test loss.}\vskip -0.5em
\label{tab:ablation}
\setlength\tabcolsep{2mm}
\begin{tabular}{l cc cc cc c}
\toprule
        & \multicolumn{2}{c}{MSFT}
        & \multicolumn{2}{c}{INTC}
        & \multicolumn{2}{c}{JPM}
        & \multicolumn{1}{c}{Avg.}\\
        \rule{0pt}{2ex}
        & \multicolumn{1}{c}{Bid} & \multicolumn{1}{c}{Ask}
        & \multicolumn{1}{c}{Bid} & \multicolumn{1}{c}{Ask}
        & \multicolumn{1}{c}{Bid} & \multicolumn{1}{c}{Ask}
        & \multicolumn{1}{c}{   } \\
\midrule
HC        &  6.79   &  7.42  &  5.47  &  6.78  &  6.22  &  3.86  & 6.09\\
ES        &  6.31   &  6.99  &  5.94  &  5.70  &  6.04  &  3.36  & 5.72\\
HC+ES     &  6.26   &  6.45  &  5.40  &  5.50  &  6.02  &  3.32  & 5.49\\
HC+ES+WS  &  6.28   &  6.15  &  5.33  &  5.64  &  6.18  &  3.38  & 5.49\\
\bottomrule
\end{tabular}
\end{table*}

Results are shown in Table~\ref{tab:ablation}.
We see that predictions from the HC alone have the highest test error, as the inputs it receives are trimmed to contain only the most relevant data for current LOB volume prediction.
ES module, by receiving complete TAQ data input and modelling market events as a stochastic process using a continuous RNN, achieves a lower error than HC module and is the dominant module that contributes to prediction accuracy. Combining both HC and ES modules, either using a predefined or adaptive weight, achieves the best performance, which suggets that the predictions from HC and ES are complementary and can be effectively combined to gain a higher accuracy prediction. The purpose of the WS module in the original LOBRM study is to facilitate the model's use in transfer learning and it does not contribute to prediction accuracy when tested on the same stock that it was trained, as is the case here.

\subsection{Superiority of Sparse Encoding for TAQ}\label{sec:encoding}
\subsubsection{LOB volume prediction} In sparse encoding, only volume information in TAQ is encoded explicitly and price information are embedded implicitly by positions of non-zero elements in one-hot vectors. In explicit encoding, information including price, volume, and trade direction are encoded directly as non-zero elements in feature vectors. Here we use the ES module in LOBRM (Decay-T) as the main model and tested the model prediction accuracy when these two different encoding methods are used. We don't choose the full model as explicitly encoded input cannot be trimmed so HC module is not used. Results are shown in Table~\ref{tab:encoding1}. We see that the model with sparse encoding achieves lower test loss error in 5 out of 6 experiments. The average test loss for the model with sparse encoding is 17.9\% lower than the model with explicit encoding.

\begin{table*}[t]
\small
\centering
\caption{LOBRM test loss comparison of explicit and sparse encodings of TAQ data.}\vskip -0.5em
\label{tab:encoding1}
\setlength\tabcolsep{2mm}
\begin{tabular}{l cc cc cc c}
\toprule
        & \multicolumn{2}{c}{MSFT}
        & \multicolumn{2}{c}{INTC}
        & \multicolumn{2}{c}{JPM}
        & \multicolumn{1}{c}{Avg.}\\
        \rule{0pt}{2ex}
        & \multicolumn{1}{c}{Bid} & \multicolumn{1}{c}{Ask}
        & \multicolumn{1}{c}{Bid} & \multicolumn{1}{c}{Ask}
        & \multicolumn{1}{c}{Bid} & \multicolumn{1}{c}{Ask}
        & \multicolumn{1}{c}{   } \\
\midrule
Explicit       &  6.87   &  6.43  &  6.80  &  9.95  &  7.60  &  4.14 &  6.97\\
Sparse         &  6.31   &  6.99  &  5.94  &  5.70  &  6.04  &  3.36 &  5.72\\
\bottomrule
\end{tabular}
\end{table*}

\subsubsection{Price trend prediction} We compare the sparse encoding method with the convolution method proposed in \cite{zhang2019deeplob} for quotes data in the task of stock price trend prediction. On the basis of explicit encoding, the convolution encoding method applies two convolution layers with filters of size $[1 \times 2]$ and stride $[1 \times 2]$ on quote data.
This structure first convolutionalize price and volume information at ask and bid sides respectively and then convolutionalize two sides' information together.
We approach a similar task presented in \cite{zhang2019deeplob} but simplify the model structure to concentrate on the prediction accuracy variation brought by different encoding methods.
We set length of time series samples as $S=50$. For sparse encoding, we set $k=5$. For convolution encoding, we use 16 $[1 \times 2]$ kernels with stride $[1 \times 2]$ followed by LeakyReLU activation. We standardize features using min-max or z-score standardization. Encoded data are passed to an MLP with ReLU activation. A GRU unit is used to receive iterative inputs and the final latent state is connected with an MLP with Softmax activation to generate a  possibility distribution over three labels $\left \{down, same, up\right \}$. We test the model on MSFT five-day dataset, using first three days for training, the fourth day for validation, and the fifth day for testing. We run a rolling average of five timesteps to alleviate label imbalance \cite{ntakaris2018benchmark}, with 29\%, 40\%, and 31\% for \textit{up}, \textit{same}, and \textit{down}. We train the model with cross entropy loss for 50 iterations and choose the model with highest validation accuracy.

Results are shown in Table~\ref{tab:encoding2}.
We can see that the sparse encoding method with two different standardization methods has superior performance in the task of price trend prediction, compared with convolution encoding either with or without price information. Thus, we draw the conclusion that the sparse encoding method for TAQ data can not only benefit the task of LOB volume prediction, but also other tasks including stock price trend prediction.

\begin{table}[t]
\small
\centering
\caption{Future price prediction: validation/test accuracy.}\vskip -0.5em
\label{tab:encoding2}
\setlength\tabcolsep{2mm}
\begin{tabular}{l cc cc cc}
\toprule
        & \multicolumn{1}{c}{Sparse}
        & \multicolumn{1}{c}{Convolution}
        & \multicolumn{1}{c}{Convolution}\\
        & \multicolumn{1}{c}{(implicit price)}
        & \multicolumn{1}{c}{(explicit price)}
        & \multicolumn{1}{c}{(no price)}\\
\midrule
Min-max       &
\multicolumn{1}{c}{59.8\% / 55.2\%}   &  \multicolumn{1}{c}{55.9\% / 54.4\%}  &  \multicolumn{1}{c}{57.2\% / 53.9\%}  \\
Z-score      &
\multicolumn{1}{c}{60.2\% / 57.4\%}   &  \multicolumn{1}{c}{57.1\% / 54.0\%}  &  \multicolumn{1}{c}{58.3\% / 55.0\%}  \\
\bottomrule
\end{tabular}
\end{table}

\subsection{Is the Model Well-Trained?}
As financial time series suffer from stochastic drift, in the sense that the distribution of data is unstable and tends to vary temporarily, large amounts of data is needed for model training. For example, the LOB used in \cite{blanchet2017unraveling} and \cite{sirignano2019universal} is of one month's and seventeen month's length.
Here, for all six sets of experiments (3 stocks $\times$ 2 sides) we use three days' LOB data for training, one day's data for validation, and one day's data for testing. We would like to test whether this amount of training data is abundant enough for out-of-sample testing.

\begin{table*}[t]
\small
\centering
\caption{LOBRM (Decay-T) test loss against training size.}\vskip -0.5em
\label{tab:historical}
\setlength\tabcolsep{2mm}
\begin{tabular}{l cc cc cc c}
\toprule
        & \multicolumn{2}{c}{MSFT}
        & \multicolumn{2}{c}{INTC}
        & \multicolumn{2}{c}{JPM}
        & \multicolumn{1}{c}{Avg.}\\
        \rule{0pt}{2ex}
        & \multicolumn{1}{c}{Bid} & \multicolumn{1}{c}{Ask}
        & \multicolumn{1}{c}{Bid} & \multicolumn{1}{c}{Ask}
        & \multicolumn{1}{c}{Bid} & \multicolumn{1}{c}{Ask}
        & \multicolumn{1}{c}{   } \\
\midrule
Day3       &  7.15   &  7.15  &  6.05  &  6.18  &  6.29  &  3.33  & 6.03\\
Day2+3     &  6.25   &  6.58  &  5.86  &  5.73  &  6.40  &  3.32  & 5.69\\
Day1+2+3   &  6.28   &  6.15  &  5.33  &  5.64  &  6.18  &  3.32  & 5.48\\
\bottomrule
\end{tabular}
\end{table*}

As the dataset we possess contains five consecutive trading days' LOB data, we leave out the fourth day's and fifth day's LOB data for validation and testing. We first use the third day's data for training and then iteratively add in the second and the first day's historical data to observe how the validation and testing loss change. Results are shown in Table~\ref{tab:historical}. We can see that there is a downward tendency in average test loss as more historical data is used for training. This is especially true for {\tt MSFT} ask, {\tt INTC} bid, and {\tt INTC} ask. This phenomenon suggests that the influence of stochastic drift may be alleviated by exposing the model to more historical samples.
Therefore, we would expect even lower out-of-sample errors if more historical data can be used for training.

\section{Conclusion}\label{sec5}
We have extended the research on the LOBRM, the first deep learning model for generating synthetic LOB data. Two major revisions were proposed: standardizing LOB data with time-weighted z-score to improve the model's generalization ability; and substituting the original ODE kernel with an exponential decay kernel (Decay-T) to improve time efficiency. Experiments were conducted on an extended LOBSTER dataset, a strict superset of the data used in the original study, with size approximately ten times larger. Using a fully chronological training and testing regime, we demonstrated that LOBRM (Decay-T) has superior performance over traditional models, and showed the efficacy of module ensembling. We further found that: (1) LOB volume prediction accuracy is negatively related to volume volatility; (2) sparse one-hot positional encoding of TAQ data  can benefit manifold tasks; and (3) there is some evidence that the influence of stochastic drift can be alleviated by increasing the number of historical samples. As a whole, this study validates the use of LOBRM in demanding application scenarios that require efficient inference and involve large amounts of data for training and predicting.

\subsubsection*{Acknowledgements}\label{ack}
Zijian Shi's PhD is supported by a China Scholarship Council (CSC)/University of Bristol joint-funded scholarship. John Cartlidge is sponsored by Refinitiv.

\bibliographystyle{splncs04}