# Generative time series models using Neural ODE in Variational Autoencoders

Marcus Lenler Garsdal, s174440  Valdemar Søgaard, s174431  Simon Moe Sørensen,
s174420  Technical University of Denmark (DTU)
02456 Deep Learning Fall 21

###### Abstract

In this paper we present an extension to the work done in [1] by implementing
Neural Ordinary Differential Equations in a Variational Autoencoder setting
for generative time series modeling. An object-oriented approach to the code
was taken to allow for easier development and research and all code used in
the paper can be found in this Github link.

The results in [1] were initially recreated and the reconstructions compared
to a baseline Long-Short Term Memory AutoEncoder. The model was then extended
with a LSTM encoder and challenged by more complex data consisting of time
series in the form of spring oscillations. The model showed promise, and was
able to reconstruct true trajectories for all complexities of data with a
smaller RMSE than the baseline model. However, it was able to capture the
dynamic behavior of the time series for known data in the decoder but was not
able to produce extrapolations following the true trajectory very well for any
of the complexities of spring data. A final experiment was carried out where
the model was also presented with 68 days of solar power production data, and
was able to reconstruct just as well as the baseline, even when very little
data is available.

Finally, the models training time was compared to the baseline. It was found
that for small amounts of data the NODE method was significantly slower at
training than the baseline, while for larger amounts of data the NODE method
would be equal or faster at training.

The paper is ended with a future work section which describes the many natural
extensions to the work presented in this paper, with examples being
investigating further the importance of input data, including extrapolation in
the baseline model or testing more specific model setups.

##  I Introduction

In the recent years, deep learning has shown great advances due to increased
computational power and availability of different frameworks. As an example,
it has shown exceptional ability in image recognition tasks using the ResNet
framework. In 2019, as an extension based on ResNet, a new framework was
introduced which uses Ordinary Differential Equations (ODEs) to represent
hidden layers in continuous time rather than discrete time, this type of
neural net is named a NeuralODE (NODE) [1]. Taking the NODE one step further,
it can be utilized in continuous time extra- and interpolation for sequential
data, by implementing the NODE as a decoder of a Variational Auto Encoder.
This paper will investigate the robustness and possibilities of this
technique.

###  I.1 From ResNet to ODE’s

ResNet utilizes ’residual blocks’ as the core part of its architecture. These
blocks are defined at a finite and discrete number of times and the output of
each block can be expressed as a function

$\mathbf{h}_{t+1}=\mathbf{h}_{t}+f(\mathbf{h}_{t},\theta_{t})$ (1)

With $t$ belonging to a finite set ranging from 0 to a fixed amount of steps
$T$, corresponding to the number of blocks in the net. Going from this
discrete setting to a continuous one can be done by looking at what happens
when the step size becomes smaller. When the step size becomes small enough,
at the limit when $t\rightarrow 0$ the change to the hidden step
$\mathbf{h}_{t+1}$ can be seen as a derivative with respect to $t$ instead
[1].

$\frac{d\mathbf{h}}{dt}=f(\mathbf{h}(t),t,\theta)$ (2)

This allows us to define the depth of the network in a continuous manner. The
function $f$ represents the hidden dynamics and is parameterized as a neural
network. To solve the derivatives given in Equation 2 the NODE uses a black
box ordinary differential equation solver, which takes as input an initial
hidden state $\mathbf{h}(0)$ to solve an initial value problem up to some
given time $T$. This way the ODE solver yields a representation of a
continuous hidden state trajectory, instead of a discrete amount of hidden
states. This also means that any specific hidden state along the hidden
trajectory can be evaluated, even with uneven step sizes, which is one of the
advantages with this approach.

###  I.2 Adjoint-state method

The training process of NODE is a bit different compared to usual Feed Neural
Networks (FFNs) which needs to store the operations of the forward pass in
order to compute the weights using the chain rule in a backpropagation
algorithm.

In the NODE framework, no intermediate operations are stored in the forward
pass (when solving the ODE of the hidden trajectory), which is what leads to
the constant memory property presented in [1]. Instead, when backpropagating
through the network $f$ to update the weights, $\theta$, a second augmented
ODE is solved, which can return the gradient of the parameters with respect to
the loss function $\frac{\partial L}{\partial\theta}$. This augmented ODE is
solved backwards in time, starting from the final time step $t_{1}$ and going
to the initial time step $t_{0}$. This is what allows the model not to store
the operations in the forward pass. The augmented ODE takes advantage of the
adjoint sensitivity method, which is a mathematical method that can
efficiently compute the gradients in an ODE. To find $\frac{\partial
L}{\partial\theta}$, a quantity called the adjoint
$\mathbf{a}(t)=\frac{\partial L}{\partial\mathbf{z}(t)}$ is defined, which
corresponds to the gradient of the loss function with respect to the hidden
state $\mathbf{z}(t)$. The derivative of the adjoint can further be expressed
as

$\frac{d\mathbf{a}(t)}{dt}=-\mathbf{a}(t)^{\top}\frac{\partial
f(\mathbf{z}(t),t,\theta)}{\partial\mathbf{z}}$ (3)

From here, the adjoint and the partial derivative of $f$ with respect to
$\theta$, can be used to find the derivative of the loss function with respect
to $\theta$.

$\frac{dL}{d\theta}=-\int_{t_{1}}^{t_{0}}\mathbf{a}(t)^{\top}\frac{\partial
f(\mathbf{z}(t),t,\theta)}{\partial\theta}dt$ (4)

This integral can be computed through a call to an ODE solver, which takes as
input the aforementioned augmented function, comprising of the hidden
trajectory $\mathbf{z}(t)$ and the adjoint state $\mathbf{a}(t)$, as well as
the initial states of the augmented function (which are defined at $t_{1}$),
namely
$s_{0}=\left[\mathbf{z}\left(t_{1}\right),\frac{\partial\bar{L}}{\partial\mathbf{z}\left(t_{1}\right)},\mathbf{0}_{|\theta|}\right]$.
It is important to note that the adjoint state depends directly on the hidden
trajectory $\mathbf{z}(t)$, but this can implicitly be computed, alongside the
adjoint state in the call to the ODE solver. If the loss depends explicitly on
previous states of the hidden trajectories (observations) it important to
adjust the adjoint state accordingly, based on the direction of the derivative
of the loss with respect to the hidden trajectory at this specific state,
$\frac{\partial L}{\partial\mathbf{z}(t_{i})}$. This is illustrated in Figure
1.

Figure 1: Reverse-mode differentation of an ODE solution[1]

###  I.3 Using Neural ODE’s as VAE

NODE’s are intrinsically well suited for temporal data and have many
advantages due to the continuous aspect of the underlying ODE. It can be used
as a generative model utilizing learned representations of a latent space. In
this setting it can be used to predict or extrapolate data, as well as
interpolate and impute missing data within a time series. This also allows the
data to be irregularly sampled when training and using the model. As a
generative model, the NODE works similarly to a VAE where a latent space
representation is learned through the training of the model. A VAE has an
encoding and decoding part, and connecting these is a ’bottleneck’ layer
containing the latent space distribution. Latent states can then be sampled
from the learned distribution and given to the decoder part of the network,
which then reconstructs (maps) the latent states into the original space as
the input data.

Our work revolves around using the NODE framework as a VAE, where it can learn
the latent space distribution of a time series. In this setting the latent
space corresponds to the initial state of the latent trajectory
$\mathbf{z_{0}}$, which an ODEsolver can take as input, in order to compute
the entire latent trajectory $z(t)$. The hidden trajectory is then evaluated
at specific time steps, and these steps are then fed to a decoder that
transforms it back into the temporal space. The workflow of the model is
sketched in Figure 2. We will try the model on different toy examples, and on
some real world solar data, and evaluate its performance against a baseline
ODE to investigate the performance and limits.

Figure 2: Structure of the neural ODE latent variable generative model[1].

##  II Methodology and Data

In this section we will discuss the methodology for how we have achieved the
results in the paper. In addition, the data used for the experiments will be
presented along with models and hyperparameters for reproducibility.

###  II.1 Experiments

We have trained the NODE model on three different data sets to test the
robustness and performance of the method across multiple data sets and use
cases. The first data set consists of clockwise and counterclockwise spirals,
the second data set consists of three subcategories of sine waves, also
referred to as springs in the paper and the third data set consists of real
world solar power curves. The data is represented as sequential data in a
matrix with size: [num data, sequence, features]. The features are (x,y)
coordinates for all data sets.

All data sets are split into a $60\%/20\%/20\%$ train/validation/test split.
Due to different complexities, the total number of data points depend on the
respective data set.

####  II.1.1 Spiral

The spiral data set is generated and consists of a total of 500 counter- and
clockwise spirals, with an even 50/50 split between each, meaning there are
250 counterclockwise spirals and 250 clockwise spirals. For each spiral (data
point) we have a total of 300 time steps where we train the models on a
subsample of 200 time steps as seen in Figure 5. Finally, to introduce some
uncertainty to the data, we introduced Gaussian noise to the subsample.

####  II.1.2 Spring

The spring data set is generated and contains three different types of
springs. A non-dampened sine wave (1), a dampened sine wave (2) and an
exponentially dampened sine wave (3) as seen in Figure 3.

Figure 3: Different spring examples included in the data: (left) example 1,
(mid) example 2 and (right) example 3.

The models are trained on 7 different combinations of the springs to
investigate how the model performs when increasing data complexity, we refer
to these combinations as mixed spring data. The combinations are:
$[1,2,3,12,13,23,123]$. For all spring data sets there are 5000 data points
with 300 time steps and trained on a subsample of 200 time steps. As an
uncertainty term, each new sampled oscillation would be adjusted by noise from
a normal distribution which would change the amplitude and shift along the
x-axis, while the frequency of the oscillation was kept constant to make it
easier for the model.

####  II.1.3 Solar power

The solar power data set comes from real life solar power data where the power
output of a solar cell is measured at a 30 minute interval throughout the day.
This yields 48 time steps. Moreover, the total number of data points is 68 and
thus not a lot of data to train on compared to the other data set.

###  II.2 Code structure

Since this paper explores a wide range of options and across multiple models
and data sets, an object-oriented approach was taken to organize code, enable
faster development and streamlined training. The code structure is elaborated
further in the GitHub README.md

###  II.3 Models

The VAE NODE models follow the structure seen in [1] and Figure 2 where time
series data is encoded into a latent state using a time variant neural net
such as a Recurrent Neural Network (RNN), or in our case a Long-Short Term
Memory (LSTM) network. The latent state is represented from
$\mathbf{z}_{t_{0}}$ to ensure that decoding happens from the beginning of the
time series. In order to enable this we encode the time series in reverse thus
ending up in $\mathbf{z}_{t_{0}}$. For the decoding process, an ODE layer is
solved using an ODE solver over a user-defined time range to compute the
latent trajectory from $\mathbf{z}(0)$ to $\mathbf{z}(t)$. This is what
enabled inter- and extrapolation in continuous time. Finally this layer is
passed to a standard FFNN with one hidden layer to map back into the data
domain.

For the baseline models, an LSTM Autoencoder is used to compare and evaluate
the reconstructions of the NODE VAE.

The baseline consists of a two-cell LSTM with a hidden layer of size 45, tanh
activation function and 2 features in the encoded dimension. This is decoded
with the model architecture but in reverse. The NODE VAE models are different
for each data set. For the spiral data set the model encodes using an RNN with
45 hidden states and tanh activation function. It has 4 latent dimensions
which are decoded using an ODE layer of a FFNN with 3 hidden layers, 30 hidden
states and an ELU activation function. Finally, a FFNN of 1 hidden layer with
30 hidden states and ReLU activation function is used to map back into the
data domain. The spring model encodes using a 1-layer LSTM with 60 hidden
states and tanh activation function. The decoding differs from the spiral
example only in 40 hidden states instead of 30. The solar power model has the
same structure as the spring model.

The baseline model used the RMSE loss and the NODE VAE models used the ELBO
loss.

It should be noted that finding the final model architectures was a highly
iterative and exploratory process as good results lead to new ideas and
improvements which in turn lead to better results and thus more ideas. Mostly
the learning rate, model width, replacing RNN with LSTM, ODE solving method
and noise in input data was configured. There is still a lot of different
configurations to try out. This will be discussed in Section V.

###  II.4 Hyperparameters

The models are trained using an ADAM optimizer with a learning rate of 0.001
for all models, and the ODE Solver utilizes the fixed-step Runge-Kutta method
(rk4). The spiral models are trained on 10k epochs, the spring models are
trained on 20k epochs and the solar power models are trained on 25k epochs. It
is also worth noting, that all training was performed using the same settings
at DTU’s High Performance Computing (HPC) clusters, to ensure that training
time etc are comparable.

##  III Results

The following section will discuss the results found from the experiments
described in the previous section. In Section III.1 we compare the
reconstructions and extrapolations of the NODE qualitatively, in Section III.2
we analyze the RMSE values and in Section III.3 we analyze the training time
per epoch across all models. Finally we visualize the latent space of one of
the trained models in Section III.4.

###  III.1 Reconstructions and extrapolation

####  III.1.1 Spiral

It can be seen from Figure 4 how the model is able to reconstruct the true
trajectory when sampled data is available, as well as extrapolate for both
$t<0$ and $t>0$ and follow the true trajectory.

Figure 4: Reconstruction and extrapolation of a single spiral sample using
NODE.

We have adjusted the extrapolation by changing the time $t$ presented to the
ODE solver, and if $t$ is set much larger or smaller than the time $t$ which
the model was trained on, then the learned trajectories begin to deviate
drastically from the ground truth. Examples of this can be seen in the
reconstruction grids in the appendix, e.g Figure 8, when both the learned
backwards and forward trajectories deviate from the true trajectory. In some
model runs the learned trajectory was experienced to explode to a completely
different order of magnitude and never follow the true trajectory. Comparing
the spiral results to [1] we find that we were not able to extrapolate as far
into the middle of the spiral as their results. However, we show that we can
extrapolate in both directions equally well, which was not visually presented
in their paper.

####  III.1.2 Springs

For visualizations of the spring data set, we present the model trained on the
Spring (1, 2, 3) data in Figure 5. It becomes clear that the model is able to
reconstruct the original trajectory when sampled data is available, however
compared to the spiral example, the extrapolation has not captured the dynamic
behavior of the oscillation for both $t<0$ and $t>0$. In the appendix in
Figure 9 it can be seen how for the simpler dynamics of example 1 the
extrapolation follows the initial direction of the oscillation, but explodes
for larger values of $t$.

Figure 5: Reconstruction and extrapolation of a single spring sample with
mixed examples 1 using NODE.

####  III.1.3 Solar power

In Figure 6 we can see that the model is able to capture the dynamics of the
solar power data in the reconstruction and generate a sinusoidal shape which
is representative of the daily power output. However, even for very small
values of $t<0$ the extrapolation can be seen to explode. It was not expected
that the model was able to extrapolate to nearby days power production, since
it was only given single-day sequences, however visually it still shows how
the extrapolation can sometimes explode to different order of magnitudes.

Figure 6: Reconstruction and extrapolation of a single solar power sample
using NODE.

###  III.2 RMSE

Table 1 shows the RMSE values for both spiral and solar power as well as
increasing complexity of spring data. For the spiral model it can be seen that
the NODE VAE model performs better than the baseline with a RMSE of 0.1284
compared to 0.1513. In the appendix Figure 11 the baseline are visualized and
it can be seen to follow the true trajectory with only slight deviations in
the beginning of the trajectory.

It is interesting that for the increasing complexity of spring data the RMSE
values are found to increase gradually from 0.0174 for Spring (1) to 0.0455
for Spring (1, 2, 3). For the baseline models the RMSE values remain more
constant around the 0.04-0.06 range across all complexities of data performing
worse than the NODE VAE. It could be investigated further if this trend
continues for the NODE VAE, and if the model will begin to struggle if
presented for even more complex data than three mixed spring examples.

For the solar data we find the RMSE for the NODE VAE and the baseline to be
very similar with respectively 0.0077 and 0.0068. Both models capture the
overall dynamic of the time-series and can reconstruct sinusoidal power curves
with the same order of magnitude as real data. None of the models are on the
other hand able to capture the intra-day variance in terms of cloud cover and
reduced power output. This might be a result of limited data availability and
could possibly be improved by including more than only 68 days of data in the
analysis.

Table 1: RMSE values for reconstruction for all model runs and baselines. Experiment | NODE | Baseline
---|---|---
Spiral | 0.1284 | 0.1513
Spring (1) | 0.0174 | 0.0647
Spring (2) | 0.0212 | 0.0582
Spring (3) | 0.0162 | 0.0471
Spring (1, 2) | 0.0261 | 0.0601
Spring (1, 3) | 0.0303 | 0.0525
Spring (2, 3) | 0.0363 | 0.0510
Spring (1, 2, 3) | 0.0455 | 0.0562
Solar Data | 0.0077 | 0.0068

###  III.3 Training time

Table 2 shows the average time per epoch for all performed experiments. The
large difference in in the magnitude of average time per epoch is due to the
changing amount of samples for training across spirals, springs and solar
data. We can see that generally the NODE is faster than the baseline for all
spring experiments, while for the spiral example they are relatively close.
For the solar data we can see quite a large difference in the average time per
epoch with the baseline being faster. This might indicate that for very little
data the NODE method does not perform so well in regards to training time,
since everything was kept constant between the spring and solar experiments
except the number of samples. It is a possibility that the NODE VAE does not
take properly advantage of its constant memory property for small amounts of
data, but this result is only an indication and something which should be
investigated further in future work.

Table 2: Average time pr. epoch for all experiments Experiment | NODE | Baseline
---|---|---
Spiral | 0.3885 s | 0.3441 s
Spring (1) | 4.9185 s | 5.7555 s
Spring (2) | 4.8733 s | 5.8242 s
Spring (3) | 4.7256 s | 5.4650 s
Spring (1, 2) | 4.7193 s | 7.5129 s
Spring (1, 3) | 4.8526 s | 5.9382 s
Spring (2, 3) | 4.9839 s | 5.9377 s
Spring (1, 2, 3) | 4.8614 s | 5.9053 s
Solar Data | 0.1467 s | 0.0533 s

###  III.4 Latent space

While evaluating the model trained on the Spring (1, 2, 3) data on its test
set, the 4-dimensional latent state $\mathbf{z}_{t_{0}}$ were sampled from the
latent space distribution $q(\mathbf{z}_{t_{0}}|\mathbf{x}_{t})$ and have been
projected in 2 dimensions using the first two principle components from a
Principle Component Analysis (PCA). The two principle components accounts for
85% of the explained variance, which should give a somewhat reliable
interpretation of the projected latent space. The 2D projections are shown in
Figure 7. The plot shows a clear distinction in the latent space between the
different type of springs used as input in the model. This is a clear
indication that the model learns to distinguish different kind of input by
sampling from different places in the distribution. Interestingly enough, we
see that the spring type 2 data, which is a damped sine wave, appears in two
separate clusters in the latent space. This is an indication that the noise we
have used when sampling the models most likely have made big enough changes to
the input data, that the model thinks spring type 2 are two separate classes
of springs. Spring type 1 has a very elongated distribution, most likely
indicating that the model has difficulties capturing the dynamics accurately.
From the RMSE of the reconstruction, it can also be seen that the spring 1
model has a higher error than the spring 3 model, which could be explained by
the model having an easier time reconstructing data that lies closer together
in the latent space. Latent space for the other models can be seen in Figure
14.

Figure 7: 2D projection of $z_{0}$ sampled from the latent space distribution
from the mixed spring 123 model.

##  IV Conclusion

To conclude on the work in this paper, we have implemented and built upon the
NODE framework presented in [1] for generative time series modeling using a
VAE framework. Several models have been trained across different complexities
of input data to test the robustness and performance of the NODE methods. In
general the NODE VAE was able to capture the dynamics of the underlying data.
We were able to reproduce the results for spirals and found that both forward
and backward extrapolation would follow the true trajectory. However, we were
not able to train the model to extrapolate for as long times as presented in
[1] (all the way to the center of a spiral).

When we subjected the NODE VAE to increasing complexity in the data, we found
a tendency that the RMSE steadily increased, while for the baseline it
remained constant. This might indicate that for even more complex data the
NODE VAE will have decreasing performance, which is an interesting hypothesis
for the method and should be investigated further in future work. The NODE VAE
was also tested in a real-life setting with little data available, and was
still able to capture general dynamics in a similar fashion to the baseline.
For all generated data sets (spiral and spring) NODE VAE outperformed the
baseline but for the real life data, the baseline performed slightly better,
in regards to reconstruction error.

Finally, when comparing the training time based on the average time per epoch,
it was found that the NODE VAE was generally faster than the baseline for
larger numbers of data. For the solar data where only 68 days of data of data
was available the baseline outperformed the NODE VAE significantly in regards
to training speed. This is most likely because of the limited data, such that
the NODE VAE does not take properly advantage of its constant memory property.

##  V Future work

Natural extensions to the work performed in this paper would include:

  1. 1.

Investigate further the importance of the input data for the NODE and how the
sampling methods affect the results.

     * •

Perform a similar analysis on samples with irregular timestamps.

     * •

Perform the analysis with Gaussian noise added to the spring data points
instead of varying the constants of the three spring equations.

     * •

Challenge the NODE model further with more complex data.

  2. 2.

Include extrapolation on the baseline to compare the forecasting capabilities
with existing methods.

  3. 3.

Extend the scope of the research by testing more specific NODE model setups.

     * •

Perform an analysis varying the ODE topology and investigating the impact on
model performance.

     * •

Test the impact of different ODE solvers.

     * •

Test the impact of different neural-net solvers (Adam, SGD).

     * •

Test how hyper-parameter tuning would impact the performance.

## References

  * Chen, Ricky T. Q. and Rubanova, Yulia and Bettencourt, Jesse and Duvenaud, David [2019] Chen, Ricky T. Q. and Rubanova, Yulia and Bettencourt, Jesse and Duvenaud, David, (2019).

## Appendix

Figure 8: Reconstruction grid for spirals visualizing multiple samples from
the NODE. Figure 9: Reconstruction grid for springs with mixed examples 1,2,3
visualizing multiple samples from the NODE. Figure 10: Reconstruction grid for
solar power visualizing multiple samples from the NODE. Figure 11:
Reconstruction grid for spirals visualizing multiple samples from the
baseline. Figure 12: Reconstruction grid for springs with mixed examples 1,2,3
visualizing multiple samples from the baseline. Figure 13: Reconstruction grid
for solar power visualizing multiple samples from the baseline. Figure 14:
Latent space of the models trained on the different combination spring data.


# Learning to Classify and Imitate Trading Agents in Continuous Double Auction
Markets

Mahmoud Mahfouz  J.P. Morgan AI Research,
Imperial College LondonLondonUnited Kingdom<EMAIL_ADDRESS>,
Tucker Balch  J.P. Morgan AI ResearchNew YorkUnited States of America
<EMAIL_ADDRESS>,  Manuela Veloso  J.P. Morgan AI ResearchNew
YorkUnited States of America<EMAIL_ADDRESS>and  Danilo Mandic
Imperial College LondonLondonUnited Kingdom<EMAIL_ADDRESS>
(2021)

###### Abstract.

Continuous double auctions such as the limit order book employed by exchanges
are widely used in practice to match buyers and sellers of a variety of
financial instruments. In this work, we develop an agent-based model for
trading in a limit order book and show (1) how opponent modelling techniques
can be applied to classify trading agent archetypes and (2) how behavioural
cloning can be used to imitate these agents in a simulated setting. We
experimentally compare a number of techniques for both tasks and evaluate
their applicability and use in real-world scenarios.

Opponent modelling, behavioural cloning, deep learning, agent-based modelling,
auction markets

††copyright: acmcopyright††journalyear: 2021††copyright:
acmcopyright††conference: 2nd ACM International Conference on AI in Finance;
November 3–5, 2021; Virtual Event, USA††booktitle: 2nd ACM International
Conference on AI in Finance (ICAIF’21), November 3–5, 2021, Virtual Event,
USA††price: 15.00††doi: 10.1145/3490354.3494386††isbn:
978-1-4503-9148-1/21/11††ccs: Computing methodologies Discrete-event
simulation††ccs: Computing methodologies Supervised learning by
classification††ccs: Applied computing Forecasting

##  1\. Introduction

The continuous double auction (CDA) is a well studied and heavily utilized
mechanism for matching buyers and sellers in a variety of markets. It is
typically used by trading exchanges in financial markets such as stocks and
commodities exchanges. A CDA is referred to as a _double_ auction as it allows
both buyers and sellers to submit orders and _continuous_ as trades are
instantaneously executed upon receipt of compatible orders (Wellman, 2011).

The financial markets can be thought of as a multi-agent system composed of a
variety of heterogeneous trading agents with different objectives and
timescales. These agents also have different utility functions (_e.g._
maximizing profit or reducing risk exposure) which are highly dependent on the
behaviour and actions of the other market participants. Exchange-traded
markets use a CDA mechanism and record the agents orders in a _Limit Order
Book_ , which we outline in section 2.1.

In this paper, we explore the feasibility of classifying and imitating the
behaviour of a representative set of trading agents in a limit order book. To
this end, we attempt to answer the following questions:

  1. (1)

Can we build a realistic agent-based model of a set of agent archetypes
trading in a limit order book?

  2. (2)

Can we learn to classify these trading agent archetypes based on the limit
order book data and agent orders?

  3. (3)

Can we imitate the behaviour of these trading agents using the limit order
book data as inputs to a model and agent order information as the labels?

Learning to classify and imitate trading agents is particularly useful for
regulators around the world who are interested in preventing the occurrence of
market manipulation. Moreover, it allows regulators to better understand the
impact of various types of participants on the market which in turn leads to
better regulations and stability in the financial markets.

##  2\. Background

###  2.1. Limit Order Books

Figure 1. Visualization of the limit order book

The limit order book (LOB) is an indicator of the supply and demand for an
exchange traded instrument at a given point in time (see figure 1). It is an
electronic record of all the outstanding buy and sell limit orders organized
by price levels. The LOB is split into two sides; the ask and bid sides
containing all the sell and buy limit orders respectively.

A matching engine is used to match the incoming buy and sell orders. This
typically follows the price/time priority rule (Preis, 2011), whereby orders
are first ranked according to their price. Multiple orders having the same
price are then ranked according to the time they were entered. If the price
and time are the same for the incoming orders, then the larger order gets
executed first. The matching engine uses the LOB to store pending orders that
could not be executed upon arrival.

Order types are distinguished between _limit_ orders and _market_ orders. A
limit order includes a price that should not be exceeded in the case of a buy,
or should not be gone below in the case of a sell. A market order indicates
that the trader is willing to accept the best price available immediately.

###  2.2. Opponent Modeling

Opponent modeling refers to a set of methods to construct models of the other
agents in the same environment with the goal of predicting various properties
of interest such as the agents’ actions, types, beliefs and goals (Albrecht
and Stone, 2018). Different approaches are studied in the literature including
policy reconstruction (Brown, 1951; Veloso, 1994), type-based reasoning
(Barrett et al., 2013; Albrecht et al., 2015), classification (Weber and
Mateas, 2009), plan recognition (Carberry, 2001), recursive reasoning (Carmel
and Markovitch, 1996), graphical models (Howard and Matheson, 2005) and group
modeling (Stone and Veloso, 1999).

Learning from data (Shoham and Leyton-Brown, 2008) is an attractive framework
under which assumptions can be relaxed and flexible models can be developed
and calibrated to real-world data. In particular, recent advances in machine
learning (Goodfellow et al., 2014; Mnih et al., 2015) allow for black-box,
task-agnostic, large-scale modeling. However, in all these settings, the
presence of multiple learning agents renders the training problem non-
stationary and often leads to unstable dynamics or undesirable final results
(Balduzzi et al., 2018; Foerster et al., 2018). Opponent modeling (Carmel and
Markovitch, 1995) addresses this problem, attempting to _decouple the sources
of non-stationarity_ by modeling separately the other learning agents from the
stationary components of the environment.

###  2.3. Behavioural Cloning

Behavioural cloning or learning from demonstrations is a supervised imitation
learning paradigm building a mapping from a historical set of observations to
actions without recovering the reward function (Osa et al., 2018). The
paradigm is particularly useful in applications whereby an agent is able to
demonstrate its behaviour rather than specify its policy directly or formulate
its reward function. Behavioural cloning has been studied and applied in
various domains including self-driving cars (Pomerleau, 1989), playing the
game Go (Silver et al., 2016), generating natural language (Wen et al., 2015)
and image captions (Karpathy and Fei-Fei, 2015). (Osa et al., 2018) provides
an excellent overview of the key methods underpinning learning from
demonstrations and current open research problems.

###  2.4. Related Work

The importance and applications of opponent modelling techniques in first-
price sealed bid iterated auctions and continuous double auction markets is
discussed in (Mahfouz et al., 2019). (Yang et al., 2012, 2015) perform
opponent modeling implicitly in a continuous double auction setup by using
Inverse Reinforcement Learning (IRL) to characterize the different trading
strategies in a limit order book. In their approach, they model the trades
placed in the order book as a Markov Decision Process (MDP) and cluster the
reward space to differentiate between high frequency, opportunistic and market
making strategies.

##  3\. Our Contributions

  1. (1)

We build an agent-based model (ABM) consisting of four different trading agent
archetypes following strategies common in practice and vary their behaviours
to ensure a wide coverage of trading strategies.

  2. (2)

We assess the realism of our market configuration and constituent agents by
comparing the statistical properties of the simulated market to those observed
in real markets.

  3. (3)

We formulate our approach to opponent modeling in a limit order book as a
supervised learning problem of identifying (classifying) agent archetypes
explicitly designed in the ABM. To this end, we designed a feed-forward neural
network classifier and compared its performance against other various
classification models.

  4. (4)

We show how a similar model to that developed for the opponent modelling task
can be used for behavioural cloning and compare the distribution of actions
between the predicted and ground truth labels.

##  4\. Agent-based Model for trading in a Limit Order Book

Agent-based modeling (ABM) is a bottom-up approach to building synthetic
multi-agent financial markets composed of various types of trading
participants. In these models, the agents are reactive to new unseen scenarios
allowing them to perform counterfactual reasoning and answer various what-if
scenarios. In our work, we utilize ABM as a method to construct a synthetic
limit order book market composed of various trading agent archetypes, which we
attempt to classify and clone their respective behaviour.

Game theory (Shoham and Leyton-Brown, 2008; Wooldridge, 2009) provides a
mathematical formalism for describing the complex behaviours that emerge from
the different agent interactions and suggests tools to develop successful
strategies. Nonetheless, strong assumptions are traditionally made about
having access to reliable models describing the system dynamics, which is
invalid even in the simplest of setups. For example, in financial markets
designing a strategy that takes into account the impact of its decisions, is
as difficult as solving the original problem itself, and heuristic approaches
(Gsell, 2008) are not adaptive hence their quality degrades as the dynamics
drift over time.

###  4.1. Simulator

We utilize and extend an open source discrete event simulation environment for
our work (Byrd et al., 2019). The simulation environment is designed for
simulating trading agents in a continuous double auction market through a
NASDAQ-like exchange agent with price-then-FIFO (first in, first out) matching
rules. It also allows for the simulation of variable electronic network
latency and agent computation delays, features that are specifically important
for high frequency trading in exchange-traded markets. The agents communicate
solely through a standardized message protocol similar to that used by NASDAQ.
The environment is controlled by a simulation kernel which manages the flow of
time and handles all the communication between the different agents in the
simulation to request latest limit order book information, place market or
limit orders, and cancel outstanding limit orders, etc.

###  4.2. Agents

Our agent population consists of a mixture of zero-intelligence (ZI) and
strategic limited-intelligence agent archetypes. The ZI agents were originally
introduced by (Gode and Sunder, 1993) to describe a family of trading agents
submitting orders in the market in a random fashion without basing the
decisions on knowledge of the current state of the limit order book.
Interestingly, it was shown in this work that the allocative efficiency of a
market arises from its structure and not the particular strategy or
intelligence of its participants.

The strategic limited-intelligence agents in our work are based on the authors
knowledge of the types of trading agents acting in real exchange-traded
markets. We refer to the zero-intelligence agents as background agents and
define three archetypes of strategic limited-intelligence agents: market
makers, market takers and directional traders which we explain in detail in
this section.

####  4.2.1. Background Agents:

Our background ZI agents consist of two agent types: (1) _Noise_ agents which
closely resemble the agents introduced in (Gode and Sunder, 1993) and (2)
_Value_ agents which base their decisions on an exogenous price time series.

The noise agents represent traders participating in the market once a day
(Kyle, 1985) by placing a market order and their participation time is sampled
from a $U$-quadratic distribution over the interval
$[t_{\text{open}},t_{\text{close}}]$, where $t_{\text{open}}$ and
$t_{\text{close}}$ are the start and end of the trading day, respectively.
This is designed to achieve higher trading activity at the beginning and at
the end of a given trading day reflecting traders’ propensity to be more
active at these times, a pattern observed in real limit order books (Dacorogna
et al., 1993).

The value agents trade based on an exogenous time series representing the
agents’ belief of the true fundamental value of the asset they are trading
(Wang and Wellman, 2017; Wah et al., 2017). In our work, we define this time
series to be the historical intra-day price of a specific stock trading in
NASDAQ.

Data: quantity $q$, simulation start time $t_{\text{sim\\_start\\_time}}$,
simulation end time $t_{\text{sim\\_end\\_time}}$, market open time
$t_{open}$, market close time $t_{close}$

NoiseAgent
_$(q,t_{\text{sim\\_start\\_time}},t_{\text{sim\\_end\\_time}},t_{open},t_{close})$_

Sample $t\sim\operatorname{UQuadratic}(t_{open},t_{close})$

if  _$t\in[t_{\text{sim\\_start\\_time}},t_{\text{sim\\_end\\_time}}]$ _ then

 _wakeup_at_ $(t)$

else

return

Sample $\epsilon\sim\operatorname{Uniform}\\{+1,-1\\}$

_place_market_order_ $(\epsilon\cdot q)$

Algorithm 1 Noise Agent trading strategy

Data: Arrival rate $\lambda_{a}$, fundamental value series
$\\{r_{t}\\}_{t>0}$, spread depth $\delta_{s}\in\mathbb{N}$, order size $q$,
inside spread probability $\xi\in[0,1]$, observation noise variance
$\sigma_{n}$

ValueAgent  _$(\lambda_{a},\\{r_{t}\\},\delta_{s},q,\xi,\sigma_{n})$_

while  __market_open_ $()$_ do

$t\leftarrow$ _current_time_ $()$

$\widehat{r}_{t}\leftarrow$ _update_fundamental_estimate_ $(r_{t},\sigma_{n})$

$b_{t},a_{t}\leftarrow$ _get_best_bid_ $()$, _get_best_ask_ $()$

$m_{t}:=\frac{1}{2}(a_{t}+b_{t})$

$s_{t}:=a_{t}-b_{t}$

Sample $X\sim\operatorname{Bernoulli}(\xi)$

if  _$X=1$_ then

$\Delta_{\text{max}}\leftarrow s_{t}\cdot\delta_{s}$

Sample $\Delta\sim\operatorname{Uniform}\\{0,1,\ldots,\Delta_{\text{max}}\\}$

else

$\Delta\leftarrow 0$

if  _ $m_{t}<\widehat{r}_{t}$ _ then

$p\leftarrow b_{t}+\Delta$

_place_limit_order_ $(q,p)$

else

$p\leftarrow a_{t}-\Delta$

_place_limit_order_ $(-q,p)$

Sample $t^{\prime}\sim\operatorname{Exponential}(\lambda_{a})$

_wakeup_at_ $(t+t^{\prime})$

Algorithm 2 Value Agent trading strategy

####  4.2.2. Market Maker Agents:

These play a central role in financial markets by providing liquidity allowing
buyers or sellers to find a match for their orders. We implement a simple
market maker in a similar fashion to that in (Wah et al., 2017; Chakraborty
and Kearns, 2011) that utilizes both the order book information and exogenous
fundamental time-series and places orders on both sides of the limit order
book with varying order sizes and participation rates.

Data: Fundamental value series $\\{r_{t}\\}_{t>0}$, Max arrival rate
$\lambda_{a}$, maximum order size $q_{\text{max}}$, maximum spread depth
$\delta_{s}\in\mathbb{N}$

MarketMakerAgent  _$(\\{r_{t}\\},\lambda_{a},q_{\text{max}},\delta_{s})$_

while  __market_open_ $()$_ do

$t\leftarrow$ _current_time_ $()$

_cancel_all_orders_ $()$

$b_{t},a_{t}\leftarrow$ _get_best_bid_ $()$, _get_best_ask_ $()$

$m_{t}:=\frac{1}{2}(a_{t}+b_{t})$

$\widehat{r}_{t}\leftarrow$ _update_fundamental_estimate_ $(r_{t})$

$p_{t}:=\frac{1}{2}(m_{t}+\widehat{r}_{t})$

sample $q_{t}\sim\operatorname{Uniform}\\{1,\ldots,q_{\text{max}}\\}$

sample $s_{t}\sim\operatorname{Uniform}\\{1,\ldots,\delta_{s}\\}$

_place_limit_order_ $(\frac{1}{2}q_{t},p_{t}-s_{t})$

_place_limit_order_ $(-\frac{1}{2}q_{t},p_{t}+s_{t})$

Sample $t^{\prime}\sim\operatorname{Exponential}(\lambda_{a})$

_wakeup_at_ $(t+t^{\prime})$

Algorithm 3 Market Maker Agent trading strategy

####  4.2.3. Market Taker Agents:

These agents place large orders over a period of time consuming liquidity from
the order book. In our experiments, we developed two types of market taker
agents adapted from common trading practices. The first TWAP agent, aims to
execute the large order at a price as close as possible to the time weighted
average price. It does this by breaking up the large order into smaller orders
of equal size and dynamically placing these orders in the market using evenly
divided time slots during a pre-defined time window. The second type of market
taker agent, VWAP agent, uses a similar approach. However, the size of the
orders placed at each point in time is weighted against the historical traded
volume of the stock. Variations of both strategies are heavily utilized in
practice by trading participants wanting to buy or sell large orders on behalf
of clients.

####  4.2.4. Directional Trader Agents:

These agents represent traders that use technical analysis to predict the
direction of the stock price. In our work, we developed two simple agents; the
first, mean-reversion agents, assume that prices revert back to a mean value
over a given time window. The second, momentum traders assume that strong
movement in prices in one direction are likely to be followed by movement in
the same direction. Both strategies implemented compare two simple moving
averages over different time periods (20 minutes and 50 minutes) and trade
accordingly. In the case of the momentum agents, the agent places a buy order
of a random size if the former moving average exceeds the latter and a sell
order otherwise.

###  4.3. Configuration and Stylized Facts

Our market configuration consists of 5100 background agents (5000 noise agents
and 100 value agents), 3 market makers, 6 market takers (3 TWAP agents and 3
VWAP agents) and 10 directional traders (5 momentum traders and 5 mean-
reversion traders). We utilize the simulator’s variable network communication
properties and use a default computational delay of 50 nanoseconds and a
pairwise latency between the agents and the exchange equivalent to the
distance between New York and Seattle. In our experiments, we use historical
mid-prices of J.P. Morgan Chase’s as the fundamental value the agents observe
and run the simulation across one week worth of NASDAQ data between 09:30 AM
and 16:00 PM which are the market open and close times of this exchange.

In order to assess the realism of our simulation, we compare a number of its’
statistical properties (referred to as _stylized facts_) to those observed in
real markets. (Bouchaud et al., 2018) provides a detailed overview of these
statistical properties covering stylized facts about asset price return
distributions, order volumes and order flow. Figures 2(a) and 2(b) show the
distribution of daily asset price returns for the simulation dataset we used
with the minutely returns displaying fat tails that become slimmer when the
period over which the returns are calculated is increased. This matches the
statistical properties observed in real market data according to (Bouchaud et
al., 2018).

(a)

(b)

Figure 2. Log returns distribution

##  5\. Classifying Trading Agent Archetypes

###  5.1. Overview

We formulate our approach to opponent modeling in a limit order book as a
supervised learning problem of identifying (classifying) agent _archetypes_
explicitly programmed in the simulator, a task which can be cast in the
auxiliary tasks framework (Jaderberg et al., 2016; Hernandez-Leal et al.,
2019) for semi-supervised reinforcement learning. Using a classification
approach for opponent modeling (Weber and Mateas, 2009) allows for predicting
properties of the other agents as opposed to predicting the future actions of
the agents, which can be useful in a range of market scenarios discussed in
our work.

###  5.2. Model

The classification model used is a simple feed-forward neural network (see
figure 3) that attempts to infer the agent archetype given the limit order
book state and the agent’s order. The state of the limit order book is given
by the z-score normalised ask and bid prices and volumes for the first 5 price
levels. The agent’s order is represented by the direction of the order (buy or
sell), its price and size. The model consists of 4 hidden layers with 256,
1024, 1024 and 1024 units respectively. In order to reduce over-fitting on the
training data, we add a drop out layer after each hidden layer. We used a
mixture of ReLU and sigmoid activation functions, an ADAM optimizer to update
the network weights based on the training data and ran the model for 200
epochs. Finally, we used an open-source hyper-parameter optimization framework
(Bergstra et al., 2013) to find the best parameters to use for our model and
dataset.

Figure 3. Neural network model architecture diagram | Market Maker | Market Taker | Directional Trader | Background Agent
---|---|---|---|---
| Precision | Recall | F-score | Precision | Recall | F-score | Precision | Recall | F-score | Precision | Recall | F-score
KNN | 0.55 | 0.69 | 0.61 | 0.47 | 0.66 | 0.55 | 0.41 | 0.29 | 0.34 | 0.41 | 0.26 | 0.32
Linear SVM | 0.73 | 0.56 | 0.64 | 0.36 | 0.24 | 0.29 | 0.41 | 0.25 | 0.31 | 0.31 | 0.61 | 0.41
RBF SVM | 0.51 | 0.80 | 0.63 | 0.61 | 0.54 | 0.57 | 0.42 | 0.35 | 0.38 | 0.44 | 0.34 | 0.38
Decision Tree | 0.99 | 0.47 | 0.64 | 0.69 | 0.70 | 0.70 | 0.55 | 0.85 | 0.67 | 0.67 | 0.64 | 0.65
Random Forest | 0.44 | 0.61 | 0.51 | 0.40 | 0.65 | 0.50 | 0.38 | 0.25 | 0.30 | 0.53 | 0.21 | 0.30
AdaBoost | 0.80 | 0.59 | 0.68 | 0.78 | 0.69 | 0.73 | 0.57 | 0.84 | 0.68 | 0.70 | 0.63 | 0.67
Naive Bayes | 0.94 | 0.47 | 0.63 | 0.18 | 0.03 | 0.06 | 0.41 | 0.92 | 0.56 | 0.57 | 0.57 | 0.57
Our model | 0.76 | 0.75 | 0.75 | 0.70 | 0.71 | 0.71 | 0.60 | 0.63 | 0.61 | 0.68 | 0.63 | 0.65
Table 1. Opponent modelling classification results

(a) KNN

(b) Linear SVM

(c) RBF SVM

(d) Decision Tree

(e) Random Forest

(f) AdaBoost

(g) Naive Bayes

(h) Our model

Figure 4. Baseline classification models confusion matrices

###  5.3. Experiments and Results

We run our ABM using fundamental values covering five trading days for one
stock and produce labelled synthetic market data that we use offline for the
opponent classification exercise. An issue arises regarding the different time
scales that each agent archetype operates on. For example, some archetypes
place substantially more orders than others, leading to the data set we
generate having high imbalance between the classes. We produce a balanced
dataset by simply down-sampling the over-represented classes.

We train our classifier on data for 3 days and evaluate its performance on the
subsequent 2 days. Figure 6 shows the confusion matrix for the four trading
archetypes in our ABM model. We compare the performance of our model against a
range of classifiers and show the precision, recall and F1-score associated
with each agent archetype in table 1.

The confusion matrix for our model shows relatively good performance when
looking at the true positives. Comparing the precision, recall and F1-scores
of our model against the baseline models shown in table 1, we see mixed
performance indicating an ensemble model is most suited. Our model performs
best in identifying the market makers and market takers and displays lower
performance identifying the directional and background traders. It is
important to note that the down-sampling of the data and the lack of long term
memory in the neural network architecture plays a significant part in this.
For example, the directional traders developed rely on past observations of
the order book to place their orders that our model and data sampling
methodology makes worse. On the other hand, we see the ensemble AdaBoost model
(Hastie et al., 2009) performing well confirming our hypothesis about the need
for an ensemble model. The AdaBoost model (Hastie et al., 2009), however, does
not show the same performance for the market makers compared to our model.

##  6\. Behavioural Cloning of Trading Agent Archetypes

###  6.1. Overview

Using the same labelled synthetic data produced by the ABM model, we formulate
our approach to behavioural cloning of the agents as a supervised regression
problem of learning the agents’ order price and size given observations of the
limit order book. This approach ties well with the paradigm taken in the
previous classification exercise and allows us to use the same dataset but
change the inputs and labels of our model.

###  6.2. Model

We use the same neural network architecture shown in figure 3 but change the
last layer of the neural network to output a vector of z-score normalized
order prices and sizes. Instead of training the model on the entire dataset,
we train it only on the observations and actions of each agent archetype. In
this model, we also used the hyper-parameter optimization framework (Bergstra
et al., 2013) to find the best parameters to use for our model and dataset for
each agent archetype.

###  6.3. Experiments and Results

We train our regression model on data specific to each agent archetype for 3
days and evaluate its performance on the subsequent 2 days. Our results
indicate the possibility of predicting the order prices and sizes with room
for improvement. Figures 5(a), 5(b), 5(c), 5(d) show the univariate
distributions of the predicted and ground truth data of both order prices and
sizes as well as a kernel density estimate for better visualisation. In almost
all cases, we see non symmetric distributions suggesting the need for a more
complex regression model.

(a) Market Makers

(b) Market Takers

(c) Directional Traders

(d) Background Agents

Figure 5. Distribution of order prices and sizes for the cloned agent
archetypes

##  7\. Conclusion

We have explored and shown the feasibility of applying opponent modelling and
behavioural cloning techniques in simulated continuous double auctions. To
this end, we built an agent-based model consisting of four different trading
agent archetypes following strategies common in practice. We assessed the
realism of our ABM against known statistical properties of real limit order
book data to ensure our simulation is realistic. We formulated our approach to
opponent modeling in a limit order book as a supervised learning problem of
identifying (classifying) agent archetypes explicitly designed in the ABM and
designed a neural network model that performs relatively well against a number
of baseline methods. Finally, we used the same data and showed the feasibility
of performing behavioural cloning by using a similar regression model to
perform this task.

The simple experiments explored and the initial results obtained are
encouraging evidence that opponent modeling and behavioural cloning methods
are potentially feasible and useful in settings of interest to the different
actors in the financial markets. Future work in this direction could explore
more unsupervised or semi-supervised learning approaches for opponent
modelling. However, the challenges anticipated following such approaches would
be around labeling any clusters obtained which would still require domain
knowledge. This suggests the need for both supervised and unsupervised
learning approaches.

Methods of opponent modeling in financial multi-agent systems can be useful in
a number of scenarios and from different perspectives, including but not
limited to:

Regulation and Safety: Market manipulation strategies represent an element of
price distortion and creation of artificial market conditions, which in turn
impedes market efficiency, transparency and fairness. In the limit order book
context, spoofing agents, for instance, submit orders only to cancel them
before their execution and then submit a buy or sell order to benefit from the
change that their initial (fake) order caused. The ability to identify these
spoofing agents and explicitly model their policies could help in the early
detection of market abuse and avoiding adverse consequences like the flash
crash that happened on May 6, 2010, causing a large and quickly reverting dip
in the Dow Jones Industrial Average stocks. Moreover, Level 4 information,
consisting of LOB labeled records, are exclusively available to exchanges and
regulators to verify that market participants are complying with rules and
processes (Paddrik et al., 2017). Having access to such a data set is an
opportunity to make a more informed and dynamic utility modeling. Quantifying
proportions of market makers, liquidity takers and buy side participants would
be valuable in detecting signs of recession or unhealthy market conditions.

Decision Making and Automation: Market participants usually only have access
to (at most) Level 3 LOB data (Paddrik et al., 2017), which is anonymized
order flow data with new, modified, and canceled limit orders placed at a
given snapshot. Although this poses a challenge to the disambiguation of limit
orders placed in the order book because it imposes a myopic approach as
opposed to using a consistent timeseries of consecutive orders placed by the
same agent, unlabeled data does not fall short of information and
opportunities. Opponent modeling can be used for the detection of iceberg
orders and other forms of hidden liquidity (Hautsch and Huang, 2012). By
studying the occurrence of such orders, a posteriori, in the LOB and the
surrounding market conditions, we can attempt to model the utility behind such
strategies. (Hautsch and Huang, 2012) studied the correlation between dark
orders and market conditions and showed that market conditions reflected by
the bid-ask spread, depth and price movements significantly impact the
aggressiveness of ‘dark’ liquidity supply. If a trader can detect such dark
orders through recognizing a pattern that iceberg order placing agents fail to
conceal or through modeling the utility of such agents and the specific cases
that make them submit such orders, they can estimate and take into account
these orders’ effect on the LOB balance and predict a larger than observed
spread.

Additionally, the Request for Quotes (RFQ) process is another setting where we
could benefit from opponent modeling. In this setting, clients send quote
requests to multiple dealers for a specific quantity of a security.
Ultimately, the dealer with the most attractive offer wins the auction. In
this single, reverse, and static auction setting, one buyer (the
auctioneer/the client) is looking to buy a single security from multiple
sellers (the dealers) who do not observe the offers submitted by the other
dealers, and only receive feedback from the auctioneer once the latter’s
decision for the winner of the auction is made. It is therefore important for
dealers receiving RFQs from different clients to have a good estimate of the
expected level at which the other dealers are willing to sell the security to
avoid mis-evaluating the offer price. This can in theory be accomplished by
having a good estimate of the fair value of the asset in question. However, it
is difficult to estimate the fair value of the asset if the asset is illiquid.
If a participant is able to infer the utilities, and by extension, the
expected offer prices of each of the other dealers in the auction, they can
narrow down their uncertainty about the fair value of the asset and avoid
mispricing it. In RFQs, apart from modeling competitors, it can be extremely
valuable to be able to model clients. In this case data is labeled, as sell
side firms (banks) are aware of their RFQ counter party (client). Performing
efficient client tiering would allow to respond to client orders in an
efficient manner, which is beneficial for both the client and the bank.

Behavioural cloning, on the other hand, is very useful in trading applications
where agents are able to demonstrate their behaviour rather than specify their
policy directly or formulate their reward. The abundance of trading agent
demonstration data stored by financial institutions makes this paradigm very
interesting to explore and presents a unique opportunity to use machine
learning for both regular and algorithmic trading applications.

###### Acknowledgements.

The authors would like to acknowledge our colleagues Thomas Spooner, Nelson
Vadori, Leo Ardon, Pranay Pasula and Rui Silva for their input and suggestions
at various key stages of the research. This paper was prepared for
informational purposes in part by the Artificial Intelligence Research group
of JPMorgan Chase & Co and its affiliates ("J.P. Morgan"), and is not a
product of the Research Department of J.P. Morgan. J.P. Morgan makes no
representation and warranty whatsoever and disclaims all liability, for the
completeness, accuracy or reliability of the information contained herein.
This document is not intended as investment research or investment advice, or
a recommendation, offer or solicitation for the purchase or sale of any
security, financial instrument, financial product or service, or to be used in
any way for evaluating the merits of participating in any transaction, and
shall not constitute a solicitation under any jurisdiction or to any person,
if such solicitation under such jurisdiction or to such person would be
unlawful.

## References

  * (1)
  * Albrecht et al. (2015) Stefano Vittorino Albrecht, Jacob William Crandall, and Subramanian Ramamoorthy. 2015.  E-HBA: Using action policies for expert advice and agent typification. In _Workshops at the Twenty-Ninth AAAI Conference on Artificial Intelligence_.
  * Albrecht and Stone (2018) Stefano V Albrecht and Peter Stone. 2018.  Autonomous agents modelling other agents: A comprehensive survey and open problems.  _Artificial Intelligence_ 258 (2018), 66–95.
  * Balduzzi et al. (2018) David Balduzzi, Sebastien Racaniere, James Martens, Jakob Foerster, Karl Tuyls, and Thore Graepel. 2018.  The mechanics of n-player differentiable games.  _arXiv preprint arXiv:1802.05642_ (2018).
  * Barrett et al. (2013) Samuel Barrett, Peter Stone, Sarit Kraus, and Avi Rosenfeld. 2013.  Teamwork with limited knowledge of teammates. In _Twenty-Seventh AAAI Conference on Artificial Intelligence_.
  * Bergstra et al. (2013) James Bergstra, Daniel Yamins, and David Cox. 2013.  Making a science of model search: Hyperparameter optimization in hundreds of dimensions for vision architectures. In _International conference on machine learning_. 115–123.
  * Bouchaud et al. (2018) Jean-Philippe Bouchaud, Julius Bonart, Jonathan Donier, and Martin Gould. 2018.  _Trades, quotes and prices: financial markets under the microscope_.  Cambridge University Press.
  * Brown (1951) George W Brown. 1951\.  Iterative solution of games by fictitious play.  _Activity analysis of production and allocation_ 13, 1 (1951), 374–376.
  * Byrd et al. (2019) David Byrd, Maria Hybinette, and Tucker Hybinette Balch. 2019\.  ABIDES: Towards High-Fidelity Market Simulation for AI Research.  _arXiv preprint arXiv:1904.12066_ (2019).
  * Carberry (2001) Sandra Carberry. 2001\.  Techniques for plan recognition.  _User Modeling and User-Adapted Interaction_ 11, 1-2 (2001), 31–48.
  * Carmel and Markovitch (1995) David Carmel and Shaul Markovitch. 1995.  Opponent modeling in multi-agent systems. In _International Joint Conference on Artificial Intelligence_. Springer, 40–52.
  * Carmel and Markovitch (1996) David Carmel and Shaul Markovitch. 1996.  Incorporating opponent models into adversary search. In _AAAI/IAAI, Vol. 1_. 120–125.
  * Chakraborty and Kearns (2011) Tanmoy Chakraborty and Michael Kearns. 2011.  Market making and mean reversion. In _Proceedings of the 12th ACM conference on Electronic commerce_. ACM, 307–314.
  * Dacorogna et al. (1993) Michael M. Dacorogna, Ulrich Alfons Müller, Robert J. Nagler, Richard B. Olsen, and Olivier V. Pictet. 1993.  A geographical model for the daily and weekly seasonal volatility in the foreign exchange market.
  * Foerster et al. (2018) Jakob Foerster, Richard Y Chen, Maruan Al-Shedivat, Shimon Whiteson, Pieter Abbeel, and Igor Mordatch. 2018.  Learning with opponent-learning awareness. In _Proceedings of the 17th International Conference on Autonomous Agents and MultiAgent Systems_. International Foundation for Autonomous Agents and Multiagent Systems, 122–130.
  * Gode and Sunder (1993) Dhananjay K Gode and Shyam Sunder. 1993.  Allocative efficiency of markets with zero-intelligence traders: Market as a partial substitute for individual rationality.  _Journal of political economy_ 101, 1 (1993), 119–137.
  * Goodfellow et al. (2014) Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio. 2014\.  Generative adversarial nets. In _Advances in neural information processing systems_. 2672–2680.
  * Gsell (2008) Markus Gsell. 2008\.  _Assessing the impact of Algorithmic Trading on markets: A simulation approach_.  Technical Report. CFS Working Paper.
  * Hastie et al. (2009) Trevor Hastie, Saharon Rosset, Ji Zhu, and Hui Zou. 2009\.  Multi-class adaboost.  _Statistics and its Interface_ 2, 3 (2009), 349–360.
  * Hautsch and Huang (2012) Nikolaus Hautsch and Ruihong Huang. 2012.  On the dark side of the market: identifying and analyzing hidden order placements.  _Available at SSRN 2004231_ (2012).
  * Hernandez-Leal et al. (2019) Pablo Hernandez-Leal, Bilal Kartal, and Matthew E Taylor. 2019\.  Agent modeling as auxiliary task for deep reinforcement learning. In _Proceedings of the AAAI Conference on Artificial Intelligence and Interactive Digital Entertainment_ , Vol. 15. 31–37.
  * Howard and Matheson (2005) Ronald A Howard and James E Matheson. 2005.  Influence diagrams.  _Decision Analysis_ 2, 3 (2005), 127–143.
  * Jaderberg et al. (2016) Max Jaderberg, Volodymyr Mnih, Wojciech Marian Czarnecki, Tom Schaul, Joel Z Leibo, David Silver, and Koray Kavukcuoglu. 2016.  Reinforcement learning with unsupervised auxiliary tasks.  _arXiv preprint arXiv:1611.05397_ (2016).
  * Karpathy and Fei-Fei (2015) Andrej Karpathy and Li Fei-Fei. 2015.  Deep visual-semantic alignments for generating image descriptions. In _Proceedings of the IEEE conference on computer vision and pattern recognition_. 3128–3137.
  * Kyle (1985) Albert S Kyle. 1985\.  Continuous auctions and insider trading.  _Econometrica: Journal of the Econometric Society_ (1985), 1315–1335.
  * Mahfouz et al. (2019) Mahmoud Mahfouz, Angelos Filos, Cyrine Chtourou, Joshua Lockhart, Samuel Assefa, Manuela Veloso, Danilo Mandic, and Tucker Balch. 2019\.  On the Importance of Opponent Modeling in Auction Markets.  _arXiv preprint arXiv:1911.12816_ (2019).
  * Mnih et al. (2015) Volodymyr Mnih, Koray Kavukcuoglu, David Silver, Andrei A Rusu, Joel Veness, Marc G Bellemare, Alex Graves, Martin Riedmiller, Andreas K Fidjeland, Georg Ostrovski, et al. 2015\.  Human-level control through deep reinforcement learning.  _Nature_ 518, 7540 (2015), 529.
  * Osa et al. (2018) Takayuki Osa, Joni Pajarinen, Gerhard Neumann, J Andrew Bagnell, Pieter Abbeel, and Jan Peters. 2018.  An algorithmic perspective on imitation learning.  _arXiv preprint arXiv:1811.06711_ (2018).
  * Paddrik et al. (2017) Mark Paddrik, Roy Hayes, William Scherer, and Peter Beling. 2017\.  Effects of limit order book information level on market stability metrics.  _Journal of Economic Interaction and Coordination_ 12, 2 (2017), 221–247.
  * Pomerleau (1989) Dean A Pomerleau. 1989\.  Alvinn: An autonomous land vehicle in a neural network. In _Advances in neural information processing systems_. 305–313.
  * Preis (2011) Tobias Preis. 2011\.  Price-time priority and pro rata matching in an order book model of financial markets.  In _Econophysics of Order-driven Markets_. Springer, 65–72.
  * Shoham and Leyton-Brown (2008) Yoav Shoham and Kevin Leyton-Brown. 2008.  _Multiagent systems: Algorithmic, game-theoretic, and logical foundations_.  Cambridge University Press.
  * Silver et al. (2016) David Silver, Aja Huang, Chris J Maddison, Arthur Guez, Laurent Sifre, George Van Den Driessche, Julian Schrittwieser, Ioannis Antonoglou, Veda Panneershelvam, Marc Lanctot, et al. 2016\.  Mastering the game of Go with deep neural networks and tree search.  _nature_ 529, 7587 (2016), 484–489.
  * Stone and Veloso (1999) Peter Stone and Manuela Veloso. 1999.  Task decomposition, dynamic role assignment, and low-bandwidth communication for real-time strategic teamwork.  _Artificial Intelligence_ 110, 2 (1999), 241–273.
  * Veloso (1994) Manuela M Veloso. 1994\.  _Planning and learning by analogical reasoning_. Vol. 886.  Springer Science & Business Media.
  * Wah et al. (2017) Elaine Wah, Mason Wright, and Michael P Wellman. 2017.  Welfare effects of market making in continuous double auctions.  _Journal of Artificial Intelligence Research_ 59 (2017), 613–650.
  * Wang and Wellman (2017) Xintong Wang and Michael P Wellman. 2017.  Spoofing the limit order book: An agent-based model. In _Proceedings of the 16th Conference on Autonomous Agents and MultiAgent Systems_. International Foundation for Autonomous Agents and Multiagent Systems, 651–659.
  * Weber and Mateas (2009) Ben G Weber and Michael Mateas. 2009.  A data mining approach to strategy prediction. In _2009 IEEE Symposium on Computational Intelligence and Games_. IEEE, 140–147.
  * Wellman (2011) Michael P Wellman. 2011\.  Trading agents.  _Synthesis Lectures on Artificial Intelligence and Machine Learning_ 5, 3 (2011), 1–107.
  * Wen et al. (2015) Tsung-Hsien Wen, Milica Gasic, Nikola Mrksic, Pei-Hao Su, David Vandyke, and Steve Young. 2015\.  Semantically conditioned lstm-based natural language generation for spoken dialogue systems.  _arXiv preprint arXiv:1508.01745_ (2015).
  * Wooldridge (2009) Michael Wooldridge. 2009\.  _An introduction to multiagent systems_.  John Wiley & Sons.
  * Yang et al. (2012) Steve Yang, Mark Paddrik, Roy Hayes, Andrew Todd, Andrei Kirilenko, Peter Beling, and William Scherer. 2012.  Behavior based learning in identifying high frequency trading strategies. In _2012 IEEE Conference on Computational Intelligence for Financial Engineering & Economics (CIFEr)_. IEEE, 1–8.
  * Yang et al. (2015) Steve Y Yang, Qifeng Qiao, Peter A Beling, William T Scherer, and Andrei A Kirilenko. 2015.  Gaussian process-based algorithmic trading strategy identification.  _Quantitative Finance_ 15, 10 (2015), 1683–1703.


# A remark of Ricci-Bourguignon harmonic soliton

Xiang-Zhi Cao School of Information Engineering, Nanjing Xiaozhuang
University, Nanjing 211171<EMAIL_ADDRESS>
###### Abstract

In this paper, we investigate the triviality of Ricci-Bourguignon harmonic
solitons. We also use the results of V-harmonic map to investigate the
property of Ricci harmonic soliton.

##  1 Introduction

Müller [21] introduced Ricci-harmonic flow, which is defined as follows: for a
closed manifold $M$, given a map $\phi$ from $M$ to some closed target
manifold $N$ :

$\frac{\partial}{\partial
t}g=-2\mathrm{Rc}+2\alpha\nabla\phi\otimes\nabla\phi,\frac{\partial}{\partial
t}\phi=\tau_{g}\phi$

where $g(t)$ is a time-dependent metric on $M,\mathrm{Rc}$ is the
corresponding Ricci curvature, $\tau_{g}\phi$ is the tension field of $\phi$
with respect to $g$ and $\alpha$ is a positive constant (possibly time
dependent).

Later, people began to study its long existence such as these works(
[6][19][18][8][17][11][25]). Fang and Zheng [13] gived the heat kernel
estimates along harmonic-Ricci flow. One can refer to the works
([20][27][7][15]) for the studies related to Ricci-harmonic flow.

Azami [4] introduced Ricci-Bourguignon harmonic flow, which is

$\displaystyle\begin{cases}\frac{\partial}{\partial t}g=-2\mathrm{Rc}-2\rho
Rg+2\alpha\nabla\phi\otimes\nabla\phi,\\\ \frac{\partial}{\partial
t}\phi=\tau_{g}\phi\end{cases}$

Next, we give a definition,

######  Definition 1.1 ( Almost Ricci-harmonic solitons).

Let $u:(M,g)\rightarrow(N,h)$ a smooth map (not necessarily harmonic map),
where $(M,g)$ and $(N,h)$ are static Riemannian manifolds.
$((M,g),(N,h),V,u,\rho,\lambda)$ is called almost Ricci-Bourguignon harmonic
solitons if

$\displaystyle\begin{cases}Rc-\rho Rg-\alpha\nabla u\otimes\nabla
u-\frac{1}{2}\mathcal{L}_{V}g=\lambda g\\\ \tau_{g}u+\langle\nabla
u,V\rangle=0,\end{cases}$

where $\alpha>0$ is a positive constant depending on $m$, $\rho$ is a real
constant and $\lambda$ is a smooth function.

When $\lambda$ is a real constant, it is called Ricci-Bourguignon harmonic
soliton . In particular, when $V=-\nabla f$, $((M,g),(N,h),f,u,\rho,\lambda)$
is called a gradient Ricci-Bourguignon harmonic soliton if it satisfies the
coupled system of elliptic partial differential equations

$\displaystyle\begin{cases}Rc-\rho Rg-\alpha\nabla u\otimes\nabla
u+\operatorname{Hess}f=\lambda g\\\ \tau_{g}u-\langle\nabla u,\nabla
f\rangle=0,\end{cases}$ (1.1)

here $f:M\rightarrow\mathbb{R}$ be a smooth tunction and where . The function
$f$ is called the potential. One can refer to [26][2][14][16][3] for the
studies on Ricci-harmonic solitons.

It is obvious that almost Ricci-Bourguignon harmonic solitons
$((M,g),(N,h),V,u,\rho,\lambda)$ is almost Ricci harmonic soliton if $\rho=0.$
Azami et al. [5] gived the condition under which complete shrinking Ricci-
harmonic Bourgainion soliton must be compact.

The gradient Ricci-harmonic soliton is said to be shrinking, steady or
expanding depending on whether $\lambda>0,\lambda=0$ or $\lambda<0$.

######  Definition 1.2.

Gradient Ricci-Bourguignon harmonic soliton is called trivial if the potential
function $f$ is constant.

It can from (1.1) that when $u$ and $f$ are constants, $(M,g)$ must be
Einstein manifold.

If $m$ is an integer with $0\leq m\leq n$ and $\alpha$ is a scalar, then the
$m$-th invariant of $\nabla^{2}f$, is denoted by $S_{m}(f)$ and is defined by
the condition (see, [22][23],pages, 461 or [24])

$\operatorname{det}\left(I+\alpha\nabla^{2}f\right)=S_{0}(f)+\alpha
S_{1}(f)+\cdots+\alpha^{m}S_{m}(f).$

Our main motivation comes from this paper [24]. On the one hand, we want to
use their methods to explore Ricci-Bourguignon harmonic soliton. On the other
hand, we know the relation between V-harmonic map and Ricci-harmonic soliton,
so we want to use V-harmonic map to study Ricci-harmonic soliton

This paper is organized as follows. In section 2, we will use the quantity
$S_{2}(f)$ to derive the triviality of Ricci-Bourguignon harmonic soliton . In
section 3, we will use the results of $V$-harmonic map to derive the property
of Ricci-harmonic solitons.

##  2 Ricci-Bourguignon harmonic soliton

In the sequel, we use the convetions and notations: $Sc=Ric-\alpha\nabla
u\otimes\nabla u-\rho Rg$, its trace is $S=(1-n\rho)R-\alpha|\nabla u|^{2}$

######  Lemma 2.1.

Let $((M^{m},g),(N,h),f,u,\rho,\lambda)$ be a gradient almost Ricci-
Bourguignon harmonic soliton, then we have the following equations:

$\displaystyle\operatorname{divSc}$ $\displaystyle=\frac{1}{2}\nabla
S-\alpha\tau_{g}(u)\nabla u,$ (2.1) $\displaystyle\langle\nabla S,\nabla
f\rangle$ $\displaystyle=2(m-1)\langle\nabla\lambda,\nabla f\rangle+2Sc(\nabla
f,\nabla f)+2\rho R|\nabla f|^{2}+(m-2)\langle\nabla R,\nabla f\rangle,$ (2.2)
$\displaystyle Sc(\nabla f,\cdot)$ $\displaystyle=\frac{1}{2}\nabla
S-(m-1)\nabla\lambda-2\rho\nabla f-(m-2)\nabla R,$ (2.3)
$\displaystyle\nabla\left(S+|\nabla f|^{2}\right)$
$\displaystyle=2(m-1)\nabla\lambda+2\lambda\nabla f+2\rho\langle\nabla
f,\cdot\rangle+(m-2)\langle\nabla R,\cdot\rangle,$ (2.4)

and

$\displaystyle\frac{1}{2}\Delta|\nabla f|^{2}=$
$\displaystyle|\nabla^{2}f|^{2}-(m-2)\langle\nabla\lambda,\nabla f\rangle-
Sc(\nabla f,\nabla f)+\alpha|\langle\nabla u,\nabla f\rangle|^{2}$ (2.5)
$\displaystyle-(m-2)\rho\langle\nabla R,\nabla f\rangle$

######  Remark 2.1.

When $((M^{m},g),(N,h),f,u,\rho,\lambda)$ is gradient Ricci-harmonic solition,
the terms about $\rho,\lambda$ and $(m-2)\nabla R$ disappear.

###### Proof.

We modify the proof of [1, Proposition 3.1] and [12, Proposition 2.1]. Using
the definition of the tensor $Sc$, we see that

$\displaystyle\frac{1}{2}\nabla R=\operatorname{div}Rc$
$\displaystyle=\operatorname{div}Sc+\alpha\operatorname{div}(\nabla
u\otimes\nabla u)+\rho g(\nabla R,\cdot)$
$\displaystyle=\operatorname{div}Sc+\alpha\tau_{g}(u)\nabla
u+\frac{\alpha}{2}\nabla|\nabla u|^{2}+\rho g(\nabla R,\cdot),$

It implies that

$\displaystyle\operatorname{divSc}$ $\displaystyle=\frac{1}{2}\nabla
S-\alpha\tau_{g}(u)\nabla u+(\frac{m}{2}-1)\rho g(\nabla R,\cdot)$

By the definition, we have

$Sc_{ij}+\nabla_{i}\nabla_{j}f=\lambda g_{ij}$

Taking trace, we get

$\displaystyle S+\Delta f=m\lambda$

Next, taking the covariant derivative me gives

$\nabla_{i}S+\nabla_{i}\nabla_{j}\nabla_{j}f=m\nabla_{i}\lambda.$

iant derivatives and using the contracted second Bianchi identity, we

$\displaystyle\nabla_{i}S=$
$\displaystyle-\nabla_{j}\nabla_{i}\nabla_{j}f+R_{il}\nabla_{l}f+n\nabla_{i}\lambda$
$\displaystyle=$ $\displaystyle-\nabla_{j}\left(-Sc_{ij}+\lambda
g_{ij}\right)+R_{il}\nabla_{l}f+n\nabla_{i}\lambda$ $\displaystyle=$
$\displaystyle\frac{1}{2}\nabla_{i}S-\alpha\tau_{g}(u)\nabla
u+(\frac{m}{2}-1)\rho g(\nabla
R,\cdot)-\nabla_{i}\lambda+R_{il}\nabla_{l}f+n\nabla_{i}\lambda$

Thus, we proved eq. 2.3

$\displaystyle\frac{1}{2}\nabla_{i}S=R_{il}\nabla_{l}f+(n-1)\nabla_{i}\lambda-\alpha\tau_{g}(u)\nabla
u+(\frac{m}{2}-1)\rho\nabla_{i}R,$ (2.6)

which is equivalent to

$\displaystyle\left(\frac{1}{2}-(n-1)\rho\right)\nabla_{i}R-\frac{1}{2}\alpha\nabla_{i}|\nabla
u|^{2}=R_{il}\nabla_{l}f+(m-1)\nabla_{i}\lambda-\alpha\tau_{g}(u)\nabla u,$

In addition, noticing that $\tau(u)=\langle\nabla u,\nabla f\rangle,$ we also
have by eq. 2.6

$\displaystyle\langle\nabla S,\nabla f\rangle$ $\displaystyle=$ $\displaystyle
2(m-1)\langle\nabla\lambda,\nabla f\rangle+2Sc(\nabla f,\nabla f)+2\rho
R|\nabla f|^{2}+(m-2)\rho\langle\nabla R,\nabla f\rangle,$

We proved eq. 2.2.

In addition,

$\displaystyle\nabla S$ $\displaystyle=2(m-1)\nabla\lambda+2Sc(\nabla
f,\cdot)+2\rho\langle\nabla f,\cdot\rangle+(m-2)\langle\nabla R,\cdot\rangle$
$\displaystyle=2(m-1)\nabla\lambda+2(\lambda g-\operatorname{Hess}f)\nabla
f+2\rho\langle\nabla f,\cdot\rangle+(m-2)\langle\nabla R,\cdot\rangle$
$\displaystyle=2(m-1)\nabla\lambda+2\lambda\nabla
f-2\operatorname{Hess}f(\nabla f)+2\rho\langle\nabla
f,\cdot\rangle+(m-2)\langle\nabla R,\cdot\rangle$
$\displaystyle=2(m-1)\nabla\lambda+2\lambda\nabla f-\nabla|\nabla
f|^{2}+2\rho\langle\nabla f,\cdot\rangle+(m-2)\langle\nabla R,\cdot\rangle,$

which implies

$\nabla\left(S+|\nabla f|^{2}\right)-2(m-1)\nabla\lambda-2\lambda\nabla
f=2\rho\langle\nabla f,\cdot\rangle+(m-2)\langle\nabla R,\cdot\rangle$

This is eq. 2.4

In the end, we prove eq. 2.5. By Bochner formula, we have

$\displaystyle\frac{1}{2}\Delta|\nabla f|^{2}$
$\displaystyle=|\operatorname{Hess}f|^{2}+\langle\nabla\Delta f,\nabla
f\rangle+Rc(\nabla f,\nabla f)$
$\displaystyle=|\operatorname{Hess}f|^{2}+m\langle\nabla\lambda,\nabla
f\rangle-\langle\nabla S,\nabla f\rangle+Rc(\nabla f,\nabla f).$

Using the second equation of the proposition yields

$\displaystyle\frac{1}{2}\Delta|\nabla f|^{2}=$
$\displaystyle|\operatorname{Hess}f|^{2}+m\langle\nabla\lambda,\nabla
f\rangle-2(m-1)\langle\nabla\lambda,\nabla f\rangle-2Sc(\nabla f,\nabla f)$
$\displaystyle\quad-2\rho R|\nabla f|^{2}-(m-2)\rho\langle\nabla R,\nabla
f\rangle+Rc(\nabla f,\nabla f)$ $\displaystyle=$
$\displaystyle|\operatorname{Hess}f|^{2}-(m-2)\langle\nabla\lambda,\nabla
f\rangle-Sc(\nabla f,\nabla f)+\alpha\nabla u\otimes\nabla u(\nabla f,\nabla
f)$ $\displaystyle-\rho R|\nabla f|^{2}-(m-2)\rho\langle\nabla R,\nabla
f\rangle,$

∎

######  Lemma 2.2 (cf.[23] ).

Let $(M,g)$ be a compact Riemannian manifold. Then the following holds:

$2\int_{M}S_{2}(f)=\int_{M}\operatorname{Ric}(\nabla f,\nabla f).$

######  Lemma 2.3 (cf. [23])).

Let $(M,g)$ be a compact Riemannian manifold. If $f$ is a smooth function on
$M$, then $S_{2}(f)$ can’t be constant unless it vanishes.

######  Theorem 2.1.

(1) If $((M^{m},g),(N,h),f,u,\rho,\lambda)$ compact gradient shrinking almost
Ricci-Bourguignon harmonic solitons with constant $S_{2}(f)$, $\rho R\leq 0$
and

$\displaystyle\int_{M}\langle\nabla\lambda,\nabla
f\rangle+\int_{M}\rho\langle\nabla R,\nabla f\rangle\leq 0$

then the soliton is trivial.

(2) If $((M^{2},g),(N,h),f,u,\rho,\lambda)$ compact gradient shrinking Ricci-
Bourguignon harmonic solitons with constant $S_{2}(f)$, then the soliton is
trivial.

###### Proof of Theorem 2.1.

As $S_{2}(f)$ is a constant, by Lemma 2.3, we infer that $S_{2}(f)$ vanishes.
Therefore, we have

$\displaystyle\int_{M}\operatorname{Ric}(\nabla f,\nabla f)=0.$ (2.7)

we obtain

$\int_{M}\langle\nabla f,\nabla
S\rangle=2(m-1)\int_{M}\langle\nabla\lambda,\nabla f\rangle+\int_{M}2\rho
R|\nabla f|^{2}+(m-2)\langle\nabla R,\nabla f\rangle.$

Now, by divergence theorem, we get

$-\int_{M}S\Delta f=2(m-1)\int_{M}\langle\nabla\lambda,\nabla
f\rangle+\int_{M}2\rho R|\nabla f|^{2}+(m-2)\langle\nabla R,\nabla f\rangle.$

Again, as $\lambda$ is a constant, we have

$\int_{M}m\lambda\Delta f+m\int_{M}\langle\nabla\lambda,\nabla f\rangle=0.$

Thus, combining , we obtain

$\int_{M}(m\lambda-S)\Delta f=(m-2)\int_{M}\langle\nabla\lambda,\nabla
f\rangle+\int_{M}2\rho R|\nabla f|^{2}+(m-2)\langle\nabla R,\nabla f\rangle.$

together with equation

$\displaystyle m\lambda-S=\Delta f,$

yields

$\int_{M}(\Delta f)^{2}=(m-2)\int_{M}\langle\nabla\lambda,\nabla
f\rangle+\int_{M}2\rho R|\nabla f|^{2}+\int_{M}(m-2)\rho\langle\nabla R,\nabla
f\rangle.$

It follows that $\Delta f=0$ in the case (1) or (2) after noticing that eq.
2.7 is

$\displaystyle\int_{M}R|\nabla f|^{2}=0.$

Since $(M,g)$ is a compact Riemannian manifold, hence $f$ is constant. This
completes the proof. ∎

Next, we give an application of the Bochner formula eq. 2.5,

######  Theorem 2.2.

(1) If $((M^{m},g),(N,h),f,u,\lambda)$ is a compact gradient Ricci-harmonic
solitons with constant $S_{2}(f)$ and $\lambda$, then this soliton is trivial.

(2) If $((M^{2},g),(N,h),f,u,\rho,\lambda)$ is a two dimensional compact
gradient almost Ricci-Bourguignon harmonic solitons with constant $S_{2}(f)$,
then this soliton is trivial.

###### Proof.

By eq. 2.5, we obtain

$\displaystyle\frac{1}{2}\Delta|\nabla f|^{2}=$
$\displaystyle\left|\nabla^{2}f\right|^{2}-(m-2)\langle\nabla\lambda,\nabla
f\rangle-Sc(\nabla f,\nabla f)+\alpha\nabla u\otimes\nabla u(\nabla f,\nabla
f)$ $\displaystyle+\rho R|\nabla f|^{2}-(m-2)\langle\nabla R,\nabla f\rangle,$

Since $S_{2}(f)$ is constant, by Lemma 2.3 we get

$\displaystyle\int_{M}\operatorname{Ric}(\nabla f,\nabla f)=0.$ (2.8)

Putting all these together and using the divergence theorem yields

$\int_{M}\bigg{(}\left|\nabla^{2}f\right|^{2}+(2-m)\langle\nabla\lambda,\nabla
f\rangle+\rho R|\nabla f|^{2}-(m-2)\rho\langle\nabla R,\nabla
f\rangle\bigg{)}=0.$

From this, we can conclude the proof easily. In the case (1), since $\lambda$
=constant, $\rho=0,$ we have $f$ is constant. In the case (2), since $m=2$,
and eq. 2.8 is just

$\displaystyle\int_{M}R|\nabla f|^{2}=0.$

∎

##  3 Ricci harmonic soliton

######  Theorem 3.1 (c.f. [9], Theorem 2.).

(1) Let $(M,g)$ be a complete noncompact Riemannian manifold with

$\operatorname{Ric}_{V}:=\operatorname{Ric}^{M}-\frac{1}{2}L_{V}g\geq-A,$

where $A\geq 0$ is a constant, $\operatorname{Ric}^{M}$ is the Ricci curvature
of $M$ and $L_{V}$ is the Lie derivative. Let $(X,h)$ be a complete Riemannian
manifold with sectional curvature bounded above by a positive constant
$\kappa$. Let $u:M\rightarrow X$ be a $V$-harmonic map such that $u(M)\subset
B_{R}(p)$, where $B_{R}(p)$ is a regular ball in $X$, i.e., disjoint from the
cut-locus of $p$ and $R<\frac{\pi}{2\sqrt{\kappa}}$. If $V$ satisfies

$\langle V,\nabla r\rangle\leq v(r),$

for some nondecreasing function $v(\cdot)$ satisfying
$\lim_{r\rightarrow+\infty}\frac{|v(r)|}{r}=0$, where $r$ denotes the distance
function on $M$ from a fixed point $\tilde{p}\in M$, then $e(u)$ is bounded by
a constant depending only on $A,\kappa$ and $R$. Furthermore, if $A=0$,
namely,

$\operatorname{Ric}^{M}\geq\frac{1}{2}L_{V}g$

then $u$ must be a constant map.

(2) Let $M^{m}$ be a complete noncompact manifold with
$\operatorname{Ric}_{f}\geq 0$ and $N$ be a complete Riemannian manifold with
nonpositive sectional curvature. If $u:M\rightarrow N$ is a f-harmonic maps
with finite weighted energy, then e(u) must be constant.

######  Theorem 3.2 (cf. [10, Theorem 12]).

Let $\left(M^{m},g\right)$ be a complete noncompact Riemannian manifold with

$\operatorname{Ric}_{V}:=\operatorname{Ric}^{M}-\frac{1}{2}L_{V}g\geq A,$

where $A\geq 0$ is a constant, $\operatorname{Ric}^{M}$ is the Ricci curvature
of $M$ and $L_{V}$ is the Lie derivative. Let $\left(N^{n},h\right)$ be a
complete Riemannian manifold with sectional curvature bounded above by a
negative constant $-\kappa^{2}(\kappa>0)$. Let $u:M\rightarrow N$ be a
$V$-harmonic map such that $u(M)\subset B_{c}$, where $B_{c}$ is a horoball
centered at $c(+\infty)$ with respect to a geodesic $c(t)$ parametrized by arc
length. Suppose that $\|V\|_{L^{\infty}(M)}<+\infty$. If
$A\geq\frac{\|V\|_{L\infty}^{2}}{m-1}$, then $u$ must be a constant map. If
$A<\frac{\|V\|_{L\infty}^{2}}{m-1}$, then $\frac{e(u)}{(B\circ u)^{2}}$ is
bounded by a constant depending only on $A,m,\kappa$ and
$\|V\|_{L^{\infty}(M)}$.

Thus, using the above Theorems, we can infer

######  Theorem 3.3.

Let u be in the same situations as Theorem 3.1, If $((M,g),(N,h),V,u,\lambda)$
is a non-compact Ricci-harmonic soliton with constant $S_{2}(f)$ and $\lambda$
, then it is non-shrinking Ricci soliton, i.e, $\lambda\leq 0$.

###### Proof.

By the assumption and Theroem 3.1, we know that $u$ is constant. So,
$\alpha$-Ricci-harmonic soliton is reduced to Ricci soliton. By contradiction,
we know that shrinking Ricci soliton is Einstein manifold (cf. [24, Theorem
1.2]) . This is a contradiction. ∎

Similarly, by the similar argument as above , we get

######  Theorem 3.4.

Let u be in the same situations as Theorem 3.2, if $((M,g),(N,h),V,u,\lambda)$
is a compact $\alpha$-Ricci-harmonic soliton with constant $S_{2}(f)$ and
$\lambda\geq\frac{\|V\|_{L\infty}^{2}}{m-1}$, then it is non-shrinking Ricci
soliton, i.e, $\lambda\leq 0$.

Recall that

######  Theorem 3.5 (cf. [28]).

Let $\left((M^{m},g),(N^{n},h),f,u,\lambda\right)$ be a shrinking or steady
gradient Ricci-Harmonic soliton satisfying with $\alpha>0$. If in addition the
sectional curvature Sect N of $N$ satisfies
$K=\sup_{N}\operatorname{Sect}^{N}<\frac{\alpha}{m}$, then $u$ is a constant
map.

Combining with Theorem 2.2 gives that

######  Corollary 3.1.

Let $\left((M^{m},g),(N^{n},h),f,u,\lambda\right)$ be a compact shrinking or
steady gradient Ricci-Harmonic soliton with $\alpha>0$, constant $S_{2}(f)$
and $\lambda$. If in addition the sectional curvature Sect N of $N$ satisfies
$K=\sup_{N}\operatorname{Sect}^{N}<\frac{\alpha}{m}$, then $(M,g)$ must be
Einstein manifold.

## References

  * [1] Abolarinwa, A. Basic structural equations for almost Ricci-harmonic solitons and applications.  Differ. Geom. Dyn. Syst. 21 (2019), 1–13.
  * [2] Abolarinwa, A., Oladejo, N. K., and Salawu, S. O. On the entropy formulas and solitons for the Ricci-harmonic flow.  Bull. Iranian Math. Soc. 45, 4 (2019), 1177–1192.
  * [3] Anselli, A. On the rigidity of harmonic-Ricci solitons.  Adv. Geom. 22, 2 (2022), 171–198.
  * [4] Azami, S. Ricci-Bourguignon flow coupled with harmonic map flow.  Internat. J. Math. 30, 10 (2019), 1950049, 17.
  * [5] Azami, S., Pirhadi, V., and Fasihi-Ramandi, G. Complete shrinking Ricci-Bourguignon harmonic solitons.  Internat. J. Math. 33, 6 (2022), Paper No. 2250046, 6.
  * [6] Băileşteanu, M. Gradient estimates for the heat equation under the Ricci-harmonic map flow.  Adv. Geom. 15, 4 (2015), 445–454.
  * [7] Băileşteanu, M., and Tran, H. Heat kernel estimates under the Ricci-harmonic map flow.  Proc. Edinb. Math. Soc. (2) 60, 4 (2017), 831–857.
  * [8] Buzano, R., and Rupflin, M. Smooth long-time existence of harmonic Ricci flow on surfaces.  J. Lond. Math. Soc. (2) 95, 1 (2017), 277–304.
  * [9] Chen, Q., Jost, J., and Qiu, H. Existence and Liouville theorems for $V$-harmonic maps from complete manifolds.  Ann. Global Anal. Geom. 42, 4 (2012), 565–584.
  * [10] Chen, Q., Jost, J., and Qiu, H. Omori-Yau maximum principles, $V$-harmonic maps and their geometric applications.  Ann. Global Anal. Geom. 46, 3 (2014), 259–279.
  * [11] Cheng, L., and Zhu, A. On the extension of the harmonic Ricci flow.  Geom. Dedicata 164 (2013), 179–185.
  * [12] Dwivedi, S. Some results on Ricci-Bourguignon solitons and almost solitons.  Canad. Math. Bull. 64, 3 (2021), 591–604.
  * [13] Fang, S., and Zheng, T. An upper bound of the heat kernel along the harmonic-Ricci flow.  Manuscripta Math. 151, 1-2 (2016), 1–18.
  * [14] Filho, W. O. C. The Poisson equation on compact Ricci solitons and Ricci-harmonic solitons.  J. Geom. 114, 2 (2023), Paper No. 14, 7.
  * [15] Guo, H., and He, T. Harnack estimates for geometric flows, applications to Ricci flow coupled with harmonic map flow.  Geom. Dedicata 169 (2014), 411–418.
  * [16] Karaca, F., and Özgür, C. Gradient Ricci-harmonic solitons on doubly warped product manifolds.  Filomat 37, 18 (2023), 5969–5977.
  * [17] Li, Y. Long time existence of Ricci-harmonic flow.  Front. Math. China 11, 5 (2016), 1313–1334.
  * [18] Li, Y. Long time existence and bounded scalar curvature in the Ricci-harmonic flow.  J. Differential Equations 265, 1 (2018), 69–97.
  * [19] Li, Y. Local curvature estimates for the Ricci-harmonic flow.  Nonlinear Anal. 222 (2022), Paper No. 112961, 53.
  * [20] Liu, X.-G., and Wang, K. A Gaussian upper bound of the conjugate heat equation along Ricci-harmonic flow.  Pacific J. Math. 287, 2 (2017), 465–484.
  * [21] Müller, R. Ricci flow coupled with harmonic map flow.  Ann. Sci. Éc. Norm. Supér. (4) 45, 1 (2012), 101–142.
  * [22] Reilly, R. C. On the Hessian of a function and the curvatures of its graph.  Michigan Math. J. 20 (1973), 373–383.
  * [23] Reilly, R. C. Applications of the Hessian operator in a Riemannian manifold.  Indiana Univ. Math. J. 26, 3 (1977), 459–472.
  * [24] Shaikh, a. a., Mandal, P., and Babu, V. a. Triviality results and conjugate radius estimation ofricci solitons.  arxive：2308.00424v1 (2023).
  * [25] Sun, J. Conditions to extend the Ricci flow coupled with harmonic map flow.  Internat. J. Math. 28, 12 (2017), 1750091, 31.
  * [26] Tadano, H. Gap theorems for Ricci-harmonic solitons.  Ann. Global Anal. Geom. 49, 2 (2016), 165–175.
  * [27] Zhu, A. Differential Harnack inequalities for the backward heat equation with potential under the harmonic-Ricci flow.  J. Math. Anal. Appl. 406, 2 (2013), 502–510.
  * [28] Zhu, M. On the relation between Ricci-harmonic solitons and Ricci solitons.  J. Math. Anal. Appl. 447, 2 (2017), 882–889.