\section{Introduction}
\noindent Following recent trends of successful AI adoption, the financial world has seen a significant surge of attempts at leveraging deep learning and reinforcement learning techniques across various application areas. Slowing down progress in this field are the particular properties of financial data: low signal-to-noise ratio ~\cite{guhr2003new}, partial observability, and irregular sampling. Furthermore, AI breakthroughs in finance often go unpublished due to monetary incentives. Additional challenges are caused by the scarcity of datasets available, which are often limited in scope, difficult to acquire or for some application areas missing altogether.

As an attempt to alleviate some of these concerns, we release both a curated dataset and a novel model for foreign exchange (FX) futures trading. We focus our attention on a particular class of FX trading methods, \textit{non-deliverable-forward} (NDF) contracts, which constitute an important open problem in finance and can serve as a challenging benchmark for supervised or reinforcement learning models. We formulate the learning problem as an optimal selection problem in which the model is tasked with selecting the end date of the forward contract (tenor) from a rich input containing past human trade patterns as well as spot rates and technical indicators. In particular, tenor selection is cast into a \textit{direct imitation learning} ~\cite{judah2012active} framework, where the model learns policy $\pi$ directly from a set of execution trajectories of a demonstration policy $\pi^*$ without receiving a reward signal from the environment. The demonstrations are derived in a greedy fashion from spot rate data and the resulting input-output tuple is utilized to perform standard supervised learning.

A key difference of our approach compared to existing FX trading algorithms lies in the type of data relied upon for learning, which includes expert tenor patterns in addition to standard technical indicators. Such patterns are extracted from a large dataset containing trades from competitive market players assumed to be informed about market state and to act rationally in order to achieve higher returns. Leveraging this additional information allows the models to differentiate between profitable and non-profitable market conditions with improved accuracy, ultimately leading to higher returns.

Fundamentally important for finance are models capable of capturing inter and intra-dependencies in highly multivariate time series. Many, if not most, of such interaction terms are nonlinear and thus challenging to analyze with standard statistical approaches. A direct consequence has been the new-found popularity of data-driven models for financial forecasting tasks, in particular recurrent neural networks (RNN) and their variants. Recurrent models, while offering an intuitive approach to time series modeling, lack an explicit module to capture inter-dependencies and perform relational reasoning ~\cite{santoro2018relational}. A different approach to time series modeling relies on temporal convolutions (TCN) as its fundamental computational block. Particularly successful in this area of research is WaveNet ~\cite{oord2016wavenet}, originally developed as a generative model for speech data. However, vanilla WaveNet and its derivative models are primarily designed to handle univariate time series and thus are ill-suited for highly multivariate financial time series. To bridge this gap, we introduce a new TCN model called \textit{WaveATTentionNet} (WATTNet) that incorporates computationally efficient dilated convolutions for temporal learning of autoregressive effects and self-attention modules to learn spatial, inter-time series interaction terms.

We summarize our main contributions as follows:
\begin{itemize}
    \item We curate, analyze, and release a new dataset containing spot rates for 64 FX currencies, along with technical indicators and hourly frequency NDF contract trade data spanning the period from $2013$ to $2019$. Several models, including classical baselines (Momentum-1, Momentum-90) and recurrent baselines (GRUs, LSTMs) are evaluated against expert NDF data.
    \item We introduce \textit{WATTNet}, a novel \textit{temporal convolution} (TCN) architecture for spatio-temporal modeling. WATTNet is designed to extend WaveNet models to settings with highly multivariate time series data.
    \item We provide two orthogonal approaches to evaluate noise stability and explain driving factors of the learned trading strategy, along with examples to highlight their efficacy.
\end{itemize}

\section{Related Work and Background}

\paragraph{Deep Learning for FX trading}
Earlier attempts at utilizing the expressivity of neural networks in \textit{forex} (FX) trading have been carried out in ~\cite{chan1995enhancing}, which predicts technical indicators via shallow fully-connected neural networks. ~\cite{yu2005designing} designs a hybrid trading system capable of providing suggestions based on qualitative expert knowledge and price forecasting data obtained from a neural network. More recently ~\cite{czekalski2015ann}, ~\cite{galeshchuk2017deep} and ~\cite{petropoulos2017stacked} leverage fully-connected neural networks, CNNs and classical autoregressive modeling techniques.
However, these approaches focus on regular forex markets and short-term predictions and rely only on spot rates and technical indicators as informative features. Incorporating additional sources of data has been explored in  ~\cite{nassirtoussi2015text}, ~\cite{vargas2017deep} ~\cite{hu2018listening}, in which additional textual information obtained from financial news articles or online discussion is included in the input features.

While the literature has no shortage of works in which reinforcement learning is applied to portfolio management ~\cite{yu2019model} and optimal trade execution, the FX markets remain comparatively unexplored. ~\cite{carapucco2018reinforcement} develops a short-term spot trading system based on reinforcement learning and obtains positive ROI in the EURUSD market. ~\cite{sornmayura2019robust} offers an analysis of deep Q-learning (DQN) performance on two FX instruments. We are not aware of any published work where deep learning or reinforcement systems are introduced to tackle FX trading in an NDF setting.

\paragraph{Spatio temporal modeling}
WaveNet ~\cite{oord2016wavenet} is an autoregressive model based on dilated temporal convolutions (TCN) in which the joint probability of the input sequence is modeled as a factorized product of conditional probabilities. SNAIL ~\cite{mishra2017simple} obtains improvements over vanilla WaveNet by adding a temporal attention layer between dilated convolutions. However, both vanilla WaveNet and SNAIL are originally designed to process univariate time series data and are thus unable to learn interaction terms between time series. ConvLSTM ~\cite{xingjian2015convolutional} introduce a convolution operation inside the LSTM cell to capture spatiotemporal information. A weakness of ConvLSTMs is given by the prior assumption of structure in the spatial domain where features closer together are prioritized by the convolution operation, as is the case for example with video data. In general applications, the time series are arbitrarily concatenated as input data and locality assumptions do not hold. Long-Short Term Network (LSTNet)~\cite{lai2018modeling} extracts local features in temporal and spatial domain with convolutions and adds a recurrent layer for longer-term dependencies. Similarly to ConvLSTM, LSTNet assumes spatial locality.
A more recent approach to spatio-temporal modeling based on Graph Neural Networks (GNNs) is Spatio-Temporal Graph Convolutional Network (STCGN) \cite{yu2018spatio} which utilizes graph convolution to carry out learning in both spatial and temporal dimensions.

\subsection{Background}
We briefly introduce the necessary background regarding different types of \textit{forex trading}.
\paragraph{Foreign Exchanges}
Trading in \textit{forex} (FX) markets is generally done via \textit{spot exchanges} or \textit{forward exchanges}, where \textit{spot rate} indicates the present expected buying rate. The \textit{spot market} can be volatile and is affected by news cycles, speculation, and underlying market dynamics. On the other hand, forward exchanges contain a long-term planning component: two parties fix a binding amount and date of exchange and the profits are calculated by comparing currency rates at the start date and fix date. The difference between start date and fix date is commonly referred to as \textit{tenor}.

\paragraph{Non-Deliverable-Forward}
An NDF operates similarly to \textit{forward exchange} contracts and exists as a replacement to forward FX trades in emerging markets. NDF markets are \textit{over-the-counter}, meaning they operate directly between involved parties without supervision, and are generally more volatile due to limited market-depth. In NDF trades the parties agree on notional amounts of primary and secondary currency (e.g. dollar USD and korean won KRW) which define the \textit{forward rate}. The currency amounts are not exchanged at the end of the contract: instead, NDF trades are cash-settled in USD, and the cash flow is computed as $R_{t,a} = (x_{t+a} - x_{t}) v_{t}$ where $x_t$ is the spot rate at time $t$, $a$ is the tenor and $v$ is the notional amount. Longer tenors are generally more profitable at the expense of a higher volatility, commonly referred to as \textit{risk-premia}. A successful trading agent thus has to find a difficult balance between risky, high return and safer, low return actions.
\section{NDF Dataset}
\paragraph{Notation}
A multivariate time series of length T and dimension M is indicated as $\{\mathbf{X}\}$. We use $\{\mathbf{x}_{i}\}$ for individual time series indexed by $i$. $x_{i,t}$ selects a scalar element of time series $i$ at time index $t$ . In particular, we indicate a slice across all time series at time $t$ with $\{x_{1,t} \dots x_{M,t}\}$.
Whenever the operations on $\{\mathbf{X}\}$ are batched, we add superscript $j$ for single samples in the batch: $\{\mathbf{X}^j\}$. With batch size $N$, the resulting tensor $\{\mathbf{X}^1\} \dots \{\mathbf{X}^N\}$ has dimensions $N \times T \times M$. We refer to tenors as $a$ and to the set of admissible tenor choices as $\mathcal{A}$.

\begin{table}
\centering
\setlength\tabcolsep{3pt}
\begin{tabular}[h!]{lr}
\toprule
NDF contracts & 7,580,814  \\
Trading hours & 35,712  \\
Trading days & 1,488  \\
Number of features per trade hour & 1,123  \\
Number of FX spot rates & 64 \\
\bottomrule
\end{tabular}
\caption{Summary of dataset statistics. Additional details included in Appendix A.}
\label{tab:datastats}
\end{table}

\paragraph{Expert benchmarks} The NDF trade records have been collected from the The Depository Trust \& Clearning Corporation (DTCC) database. These records contain start and end dates of each NDF contract, along with currency amounts. For each trading day $t$ and admissible tenor $a \in \mathcal{A}$ we obtain trading volumes $v_{t,a}$. We refer to \textit{Expert} as a trading agent that chooses tenors corresponding to maximum volumes $a_t = \argmax_{a \in \mathcal{A}}{v_{t,a}}$. In addition to \textit{Expert} we obtain a fictitious agent based on NDF records which is assumed to have partial future knowledge of the market dynamics, which we refer to as \textit{Expert oracle}. \textit{Expert oracle} is a filtered version of
\textit{Expert}: at each trading day $t$ it selects the shortest tenor with positive return:
\begin{equation}
    a_t = \argmin_{a \in \mathcal{A}}{\{ a | (x_{t+a} - x_t) > 0}\}
\end{equation}
In particular, \textit{Expert oracle} is designed to select shortest tenors to avoid a perfect accuracy exploitation of risk-premia which would set an unrealistic benchmark for any model.
\textit{Expert} and \textit{Expert oracle} are used as human trader benchmarks.
\paragraph{Multivariate input data}
In order to learn how to effectively choose NDF tenors, the models have access to time series that can be broadly categorized in three groups:
\textit{FX spot rates}, \textit{technical indicators}, and \textit{NDF tenor volumes}. Daily FX spot rates serve as contextual market information and provide a frame of reference that aids the model in identifying profitable states. We include spot rates of 64 major and minor FX currency pairs, with a complete list provided in Appendix A.
Raw financial market data is often augmented with hand-crafted features to help combat noise and non-stationarity ~\cite{ntakaris2019feature}. To this end, we choose the following technical indicators and include details in Appendix A:
\begin{itemize}
    \item Simple moving average (SMA)
    \item Exponential moving average (EMA)
    \item Moving Average Convergence Divergence (MACD)
    \item Rolling standard deviation (RSD)
    \item Bollinger Bands (BB)
    \item ARIMA 1-day spot rate forecast
\end{itemize}
The last category of input features, \textit{NDF tenor volumes}, is obtained from DTCC NDF records. For each NDF pair under consideration and each admissible tenor $a \in \mathcal{A}$, we generate a time series of volumes $v_{t,a}$ which includes a summation over all NDF records at a specific day $t$. In particular, given a choice of maximum tenor of 90 days, each NDF pair contributes with a 90-dimensional multivariate volume time series to the input, which further emphasizes the need for a model capable of processing and aggregating information across highly multivariate time series. The code for downloading and preprocessing the data will be released after publication.

\begin{figure*}[!h]
   \centering
    \includegraphics[scale=0.50]{watt_new_cmyk.pdf}
    \caption{WATTNet overview: dilated TCNs are independently applied to each univariate input time series. A single dot-product attention head subsequently aggregates information across slices $\{z_{1,t} \dots z_{M,t}\}$ and the result $\{\tilde{\textbf{Z}}\}$ is passed to the next WATTBlock or used directly for the task at hand.}
    \label{fig:watt}
\end{figure*}

\section{Model}
\paragraph{Spatio-temporal modeling with WATTNet}
\textit{WaveATTentionNet} (WATTNet) is a novel model designed for highly multivariate time series inputs. WATTNet includes \textit{temporal modules}, tasked with independently aggregating information across time steps of univariate time series $\{\mathbf{x}\} \in \{\mathbf{X}\}$ and \textit{spatial modules} which aggregate features across slices of all time series at a specific time $t$ $\{x_{1,t}, \dots x_{M,t}\}$. Temporal and spatial modules are alternated and allow for learning a hierarchical spatio-temporal representation. An overview of the model is given in Figure~\ref{fig:watt}.

\paragraph{Temporal learning}
Temporal learning is achieved by applying temporal dilated convolutions (TCN) to univariate time series $\{\mathbf{x}\} \in \{\mathbf{X}\}$. In particular, given a convolution with kernel size $k$ and dilation coefficient $d$, we compute the output at time $t$ of a dilated convolution of $\{\mathbf{x}\}$ as:
\begin{equation}
    z_t = \sum_{i=1}^{k} w_i * x_{t-i*d}
    \label{eq:dil}
\end{equation}
where $w_i$ is the $i^{\text{th}}$ weight of the convolutional kernel. Each univariate time series has access to its own set of convolutional weights $\mathbf{w}$ as temporal convolution operations are carried on independently. We note that independence between convolutions is necessary to provide the model with enough flexibility to treat time series with different characteristics. The outputs of the TCN operation are then concatenated as to form a multivariate latent time series $\{\mathbf{Z}\}$. In particular, WATTNet includes gated convolutions, a standard architectural component for sequential data. Two dilated TCNs are applied to $\{\mathbf{X}\}$ and the results $\{\mathbf{Z}\}_{\alpha}$, $\{\mathbf{Z}\}_{\beta}$ are passed to non-linear activation functions and then multiplied element-wise:
\begin{equation}
    \{\mathbf{Z}\} = \sigma(\{\mathbf{Z}\}_{\alpha}) \odot  \text{tanh}(\{\mathbf{Z}\}_{\beta})
    \label{eq:dil}
\end{equation}
where $\sigma$ indicates a sigmoid activation. The output $\{\mathbf{Z}\}$ is then fed into a spatial learning module after which the process repeats for a number of times depending on WATTNet's layer depth.

\paragraph{Spatial learning}
A single-head scaled-dot product attention mechanism ~\cite{vaswani2017attention} is placed between dilated TCN layers and allows the model to exchange information across different input time series at a specific time slice. We compute key $\mathbf{K}$, query $\mathbf{Q}$ and value $\mathbf{V}$ by considering a slice $\{z_{1,t} \dots z_{M,t}\} \in \{\mathbf{Z}\}$ of latent time series at time $t$ as the input of learnable linear transformation of the type $\psi(\{z_{1,t} \dots z_{M,t}\})$ with weights $\mathbf{W}_k$, $\mathbf{W}_q$, $\mathbf{W}_v$. The resulting matrices $\mathbf{K}$, $\mathbf{Q}$, and $\mathbf{V}$ are then used in the standard scaled-dot product attention to return M weighted averages of values $\mathbf{V}$:
\begin{equation}
\{\tilde{z}_{1,t} \dots \tilde{z}_{M,t}\} = \text{softmax}(\frac{\mathbf{Q}\mathbf{K}^T}{\sqrt{d_k}}\mathbf{V})
\end{equation}
where $d_k$ is a scaling factor given by the second dimension of $\mathbf{K}$.
The process is repeated for latent feature slices $\{z_{1,t} \dots z_{M,t}\}$, $t = 1, \dots T$ and the results are concatenated into $\{\tilde{\mathbf{Z}}\}$, a spatio-temporal latent representation of input data $\{\mathbf{X}\}$. Weights $\mathbf{W}_q, \mathbf{W}_k, \mathbf{W}_v$ are shared across the entire sequence length $T$, allowing the attention head to capture time-invariant features that incorporate information from multiple time series. Output $\{\tilde{\mathbf{Z}}\}$ can be used directly for different tasks to perform decision making conditioned on multivariate time series data; alternatively, if the task at hand benefits from deeper models, $\{\tilde{\mathbf{Z}}\}$ can instead be passed to the following TCN layer to perform additional cycles of temporal and spatial learning.

\paragraph{Hierarchical representation}
A single temporal and spatial module constitute a full WATTNet layer of computation and is referred to as \textit{WATTBlock}. WATTBlocks can be stacked, in which case output $\{\tilde{\mathbf{Z}}\}$ becomes a hierarchical spatio-temporal representation of $\{\mathbf{X}\}$. As is the case with other TCN-based models, the dilation coefficient is doubled each temporal module as to provide an increasing receptive field which allows for a computationally inexpensive way to model long sequences. An additional benefit of the gradual dilation increase is the slow introduction of interaction terms between time series which include less lagged values for early WATTBlocks and more for later ones. At layer $i$, the dilated TCN for scalar output $z_t$ has a receptive field of $2^{i}k$, with $k$ being the size of the convolutional kernel. During spatial learning, the information flow across a slice of latent TCN output $\{\mathbf{Z}\}$ at time $t$ $\{z_{1,t}, \dots z_{M,t}\}$ is thus limited to $2^{i}k$ lagged values of the raw $\{\mathbf{X}\}$, given by $\{x_{1,[t-2^{i}k:t]}, \dots x_{M,[t-2^{i}k:t]}\}$. We observe that gradually increasing the size of this interaction window is key in learning a hierarchical representation of the data that strongly intertwines spatial and temporal causal effects.

\paragraph{Computational requirements}
WATTNet can be used as a general lightweight tool for spatio-temporal modeling. The temporal modules are fully parallelizable due to complete independence of inputs and convolutional weights across different univariate time series.
To leverage an already existing fast CUDA implementation of parallel TCNs in PyTorch, we utilize \textit{grouped-convolutions}. In particular, the $M$ dimension of input time series $\{\mathbf{X}\}$ becomes the channel dimension for the TCN, and different convolutional kernels are applied to each input channel to obtain the corresponding output.
On the spatial learning front, the attention modules have a computational cost of $O(TM^2)$, which is comparable to the standard quadratic attention cost of $O(M^2)$ when $T << M$.

\section{NDF Tenor Selection}
Selecting a profitable tenor is challenging since it burdens the model with a choice of short tenors with smaller returns or long, risky tenors with a potentially greater return. One approach to training a tenor selection model is performing imitation learning on \textit{Expert} or \textit{Expert oracle} labels. Both have advantages and disadvantages; training on \textit{Expert} allows for daily online training and thus reduces the need for the model to extrapolate to periods further into the future. This aspect can be particularly beneficial for turbulent markets that display frequent regime switches. \textit{Expert oracle} labels, on the other hand, require information from up to $N$ days in the future, with $N$ being the maximum allowed tenor, since positive return filtering can only be performed by leveraging spot rate data. \textit{Expert oracle} labels can be advantageous since they teach the model to be \textit{risk-averse}; however, both approaches potentially include unwanted human-bias in the learned strategy. We propose an alternative approach in which WATTNet learns from optimal greedy tenor labels obtained directly from market data. Given the spot rate value for target FX pair $\{\mathbf{y}\}$, we extract the optimal tenor label $\bar{a}_t$ at time $t$ as:
\begin{equation}
    \bar{a}_t = \argmax_{a\in\mathcal{A}}{(y_{t+a} - y_{t})}
    \label{e}
\end{equation}
\textit{Policy divergence} ~\cite{ross2011reduction} is a performance degrading issue often present in imitation learning where the agent accumulates small expert imitation errors along the rollout and ends up in unexplored regions of the state space. To sidestep this issue, we base our strategy on a conditional independence assumption between tenor $a_t$ and state dynamics $s_{t+1}$: $p(s_{t+1}|s_{t}) = p(s_{t+1}|s_{t},a_{t})$.
The training procedure is then carried out as follows. Input time series $\{\mathbf{X}\}$ is split into overlapping length $T$ slices which are then shuffled into a training dataset. At each trading day $t$, the model $\phi$ is trained via standard gradient-descent methods to minimize the cross-entropy loss of outputs $\phi(\{\mathbf{X}_{[t-T:t]}\})$ and tenor labels $\bar{a}_t$.

\begin{table*}[h!]
\small
\centering
\setlength\tabcolsep{3pt}
\begin{tabular}{lrrrrrrrrr}
\toprule
 && USDCNY &&& USDKRW &&& USDIDR & \\
\cmidrule(r){1-10}
Model & ROI & opt.acc & nn.acc & ROI & opt.acc. & nn.acc. & ROI & opt.acc. & nn.acc. \\
\midrule
Optimal & 759.8 & 100 & 100 & 844.2 & 100 & 100 & 1260.0 & 100 & 100 \\
Expert (oracle) & 77.6 & 21.6 & 100 & 139.4 & 13.7 & 100 & 152 & 3.6 & 100 \\
\midrule
Expert & 0.0 & 1.4 & 47.8 & 12.7 & 1.9 & 49.9 & 230 & 0.4 & 67.0 \\
Momentum-1 & 14.6 & 1.4 & 48.3 & 10.7 & 0.9 & 49.7 & 201 & 0.7 & 66.5 \\
Momentum-90 & 4.9 & 6.7& 54.1 & 31.9 & 2.6 & 56.2 & $\textbf{338}$ & $\textbf{1.9}$ &69.2 \\
GRU-I  & $26.7.\pm48.5$ & $4.9\pm0.7$ & $54.7\pm1.0$ & $-98.7\pm35.4$ & $1.5\pm0.5$ & $52.1\pm1.4$ & $83.5\pm33.1$ & $0.6\pm0.1$ & $62.9\pm2.6$  \\
LSTM-I & $74.3\pm37.3$ & $3.7\pm0.9$ & $58.7\pm2.2$ & $74.6\pm30.0$ & $2.6\pm0.6$ & $56.0\pm2.5$ & $146.4\pm40.4$ & $1.1\pm0.3$ & $66.5\pm1.1$ \\
WATTNet & $\textbf{219.1}\pm25.5$ & $\textbf{6.7}\pm0.7$ & $\textbf{59.5}\pm1.6$ & $\textbf{142.4}\pm16.9$ & $\textbf{2.7}\pm0.2$ & $\textbf{59.5}\pm1.0$ & $280.2\pm39.9$ & $1.3\pm0.3$ & $\textbf{69.5}\pm0.9$ \\
\bottomrule
\end{tabular}
\caption{Test results in percentages (average and standard error). Best performance is indicated in bold.}
\label{tab:allresone}
\end{table*}

\begin{table}
\small
\centering
\setlength\tabcolsep{3pt}
\begin{tabular}[h!]{lrrrr}
\toprule
Market & $\mu_{\text{train}}$ & $\sigma_{\text{train}}$ & $\mu_{\text{val}}$ & $\sigma_{\text{val}}$\\
\midrule
USDCNY & 11.33 & 231.63 & -2.96 & 232.88 \\
USDIDR & 22.04 & 506.17 & 3.48 & 335.75 \\
USDINR & 2.88 & 297.40 & 8.19 & 286.79 \\
USDKRW & 0.50 & 471.25 & -12.75 & 387.49 \\
USDPHP & 10.44 & 271.37 & -0.79 & 240.37 \\
USDTWD & -3.73 & 307.28 & -2.09 & 235.29 \\
\bottomrule
\end{tabular}
\caption{Statistics of 1-step (daily) percent returns. Units reported are scaled up by $1e^{5}$ for clarity.}
\label{tab:prmet}
\end{table}
\begin{figure}[t]
    \centering
    \includegraphics[scale=0.19]{actions_watt_cmyk.pdf}
    \caption{Tenor actions across USDKRW test set. Background gradient is used
    to indicate the ROI of different tenors. The gradient is
    slanted since the raw return at day $t$ with a tenor of 90 days
    is the same as the return from trading day $t + i$ with tenor of
    $90 - i$ days. WATTNet learns to correctly exploit long periods of positive trend and places less actions than LSTM-I in negative (red) return regions.}
    \label{fig:actions}
\end{figure}
\section{Experimental Results}
\paragraph{NDF Markets}
The experimental evaluation covers the following 6 major NDF markets: \textit{Chinese Yuan} (USDCNY), \textit{Indonesian Rupiah} (USDIDR), \textit{Indian Rupee} (USDINR), \textit{Korean Won} (USDKRW), \textit{Philippine Peso} (USDPHP), \textit{Taiwan Dollar} (USDTWD). We elaborate on the results for USDCNY, USDKRW, USDIDR and include discussion of the remaining 3 markets as supplementary material (Appendix C). The selection has been carried out to test the proposed method on markets with different characteristics as shown in Table \ref{tab:prmet}.

\paragraph{Recurrent and classical baselines}
Two-layer stacked GRU and LSTMs are used as additional baselines for tenor and referred to as GRU-I and LSTM-I. Both recurrent baselines and WATTnet are connected to a fully-connected head which takes as input the latent representation of $\{\mathbf{X}\}$ produced by the model. A probability distribution over tenor actions is then obtained via softmax. Additionally, we include the following classical trading baselines:
\begin{itemize}
    \item \textit{Momentum-1}: 1-day lag of expert tenor actions. Effective in markets where monotonic behavior in the spot rate is frequent and the monotonic sequences span several trading periods.
    \item \textit{Momentum-90}: best performing tenor from 90 days prior. Effective in markets with trends whose duration is longer compared to the maximum tenor.
\end{itemize}

\paragraph{Training setup} The models are implemented in PyTorch and trained using Adam \cite{kingma2014adam} and a learning rate cosine decay schedule from $6e^{-4}$ down to $3e^{-4}$. To avoid overfitting uninformative noisy patterns in stale data input sequence length is set to $30$ days. In addition, to enable a fair comparison and avoid additional overfitting we employ an early stopping scheme based on training loss that is motivated by different convergence times of different models. We use a static testing approach with a long period of $446$ out-of-sample trading days to test stability of the learned trading strategy under turbulent market conditions and a wider distributional shift between in-sample and out-of-sample data. PyTorch code for models and training procedure is included in the supplementary material.

\paragraph{Metrics}
The following metrics are used to benchmark the performance of trading models and baselines:
\begin{itemize}
    \item \textit{Return on investment} (ROI): given a tenor action $a$ at time $t$ and spot rate value $x_t$, the percent ROI is calculated as $\text{ROI}_t = 100 \left( \frac{x_{t+a} - x_t}{x_t}\right)$
    \item \textit{Optimal accuracy}: standard supervised learning accuracy of model outputs $\hat{y}_t$ versus \textit{optimal} tenor labels $y_t$.
    \item \textit{Non-negative return accuracy}: accuracy of model outputs $\hat{y_t}$ compared to tenor actions with positive or zero return. At time $t$, there are generally multiple tenor actions with non-negative return, thus rendering \textit{non-negative accuracy} a less strict metric compared to \textit{optimal accuracy}. It should be noted that it is possible for a high ROI trading model to show poor optimal accuracy but competitive positive return accuracy since non-negative accuracy also captures positive ROI strategies that differ from optimal tenor labels.
\end{itemize}

\paragraph{Discussion of results}
We characterize the 6 NDF markets under evaluation based on mean $\mu$ and standard deviation $\sigma$ of their 1-day returns and performance of classical baselines. Mean-variance statistics given in Table \ref{tab:prmet} show easier markets with similar volatility in training and test periods (e.g. USDINR) as well as markets that appear more erratic and thus challenging to trade profitably in. From Table \ref{tab:allresone}, we determine USDCNY and USDKRW to be challenging for \textit{Momentum} and recurrent baselines, in addition to being barely profitable for expert traders. GRU-I is unable to consistently get positive ROI and we suspect this is caused by its insufficient ability to exchange information between time series. LSTM-I, on the other hand, fares better by leveraging its memory module to perform rudimentary spatial learning ~\cite{santoro2018relational}. USDIDR appears to be highly profitable for \textit{Momentum-90}, a phenomenon highlighting a longer positive trend and a higher density of positive return trades in the test data.
\begin{figure*}[h!]
   \centering
    \includegraphics[scale=0.25]{input_grads_cmyk.pdf}
    \caption{Highest impact spot rate features behind 90 day tenor actions in USDCNY testing period. Background color shows $20$ day rolling standard deviation.}
    \label{inputgrad}
\end{figure*}
\begin{figure}[h!]
   \centering
    \includegraphics[scale=0.34]{umap_paper_two.pdf}
    \caption{UMAP embedding of model latents. The points are labeled according to the final tenor output of the model.}
    \label{fig:umap}
\end{figure}

\begin{itemize}
    \item \textbf{USDCNY:} The results in Table \ref{tab:allresone} show that all classical baselines, including expert actions, perform poorly in terms of ROI and optimal tenor accuracy. LSTM and WATTNet, on the other hand, are able to generalize better, with WATTNet surpassing ROI of oracle trades.
    \item \textbf{USDKRW:} The USDKRW train-test split presents the most challenging trading period across all 6 NDF markets considered in this work. The return statistics in Table \ref{tab:prmet} show a decrease in mean return as well as a significant change in volatility. Traditional baselines perform poorly whereas WATTNet surpasses \textit{Expert oracle} ROI. Figure~\ref{fig:actions} highlights the ability of WATTNet to adjust its tenor actions depending on trading spot rate trends.
    \item \textbf{USDIDR:} USDIDR has been chosen to evaluate the baselines under a profitable trading period. All classical baselines trade with positive ROI, even surpassing oracle ROI due to their propensity to trade at longer tenors. WATTNet achieves performance competitive with \textit{Momentum} baselines, showing that it is capable of fully exploiting trading periods with long positive trends and a wide shift in volatility between training and test data.
\end{itemize}
\subsection{Explainability}
Model explainability is particularly important in application areas where the models are tasked with critical decision making, as is the case for algorithmic trading. Understanding driving factors behind a trading decision is necessary to properly assess the risks involved. We tackle this issue by proposing two orthogonal approaches for evaluating the tenor selection strategy in terms of noise stability and driving factors.
\paragraph{Feature importance by input gradients}
To pinpoint the driving factors of trades at different tenors we propose sorting the features by their gradient magnitude. In the case of tenor selection, each input feature carries a specific meaning which can be leveraged by domain experts to confirm whether the model outputs actions consistent with market dynamics.

Given a batch of $N$ multivariate input sequences $\{\mathbf{X}^1\} \dots \{\mathbf{X}^N\}$ with tenor labels equal to $a$, we compute the cross-entropy loss of model $\mathcal{L}$ and derive the empirical expectation for the absolute value of time series $\{\mathbf{x}_j\}$ input gradient as follows:
\begin{equation}
 \mathcal{G}_j = \frac{1}{T} \left|\sum_{t=1}^{T}\frac{\partial \mathcal{L}}{\partial x_{jt}}\right|
\end{equation}
To illustrate the effectiveness of this approach we select $6$ spot rate features with highest absolute gradient values for tenor actions of $90$ days in the USDCNY test data: EURSGD, GBPAUD, USDCHF, EURDKK, USDSGD, AUDNZD (Figure~\ref{inputgrad}). The Pearson's correlation coefficient $\rho$ between USDCNY and each of the above-listed features is computed with training and test sets. In the background, 20 day rolling standard deviation highlights regions of low and high volatility. Input sequences $\{\mathbf{x}\} \in \{\mathbf{X}\}$ which are mapped by the model to $90$ tenor actions are colored in red, whereas the actions themselves are indicated as black dots. The model learns to trade on long tenors when currencies that are positively correlated with USDCNY, such as USDSGD, undergo periods of growth. The degree to which such trends affect the model is directly reflected in $\rho$: USDCHF, still positively correlated with USDCNY, shows a less decisive positive trend, with additional ups and downs. Moreover, the model learns to favor trading periods with low volatility.
This analysis can be extended by domain experts to additional input features, such as technical indicators or past tenor actions, and can boost confidence in the decisions made by the model.

\paragraph{Latent space representation}
Desired properties of the learned trading strategy are \textit{input coherence} and \textit{stability}.
Input coherence is achieved by a model that outputs similar tenors for similar states. Stability, on the other hand, is concerned with how much noise perturbation is required to cause a tenor switch from a certain state.
We perform a visual inspection of these properties via a \textit{uniform manifold approximation and projection} (UMAP) which excels at capturing both local and global structure of the high-dimensional data ~\cite{mcinnes2018umap}. For each model, latent vectors of their last layer are embedded into two-dimensional space. UMAP outputs compact clusters of labels for input coherent models and more volumetric clusters for stable models. From Figure~\ref{fig:umap} we observe that WATTNet learns a coherent latent representation that clusters low and high tenors correctly and covers a larger volume of the embedding space. Instability of GRU-I and LSTM-I observed in the results of Table \ref{tab:allresone} can in part be explained by noticing that their learned representation lies on thin lower-dimensional manifolds with mixed tenor labels. As a result, small noise perturbations can cause wide jumps in tenor actions, potentially causing a drop in performance.

\section{Conclusion}
We introduced a challenging imitation learning problem, tenor selection for \textit{non-deliverable-forward} (NDF) contracts. With the goal of promoting further research in this direction, we constructed and released a comprehensive NDF dataset and designed \textit{WaveATTentioNet} (WATTNet), a novel model for spatio-temporal data which outperforms expert benchmarks and traditional baselines across several NDF markets. Finally, we employed two explainability techniques to determine driving factors and noise stability of the learned tenor strategy. Future work includes defining and solving \textit{order sizing} of an NDF contract, as well as augmenting the tenor selection model with a reinforcement learning agent head. Reminiscent of world models ~\cite{ha2018recurrent}, such an approach would sacrifice explainability for additional flexibility; the agent could for example be equipped with a memory module consisting of a recurrent network to keep track of its previous trades. A different direction would be exploring the use of multi-head attention instead of single-head as an effective modification to WATTNet in the case of long input sequences or multiple relations between features.

{\small
\selectfont
\bibliographystyle{aaai}


\section{Introduction}

Quantum chromodynamics (QCD) predicts that strongly interacting matter undergoes a phase
transition to quark-gluon plasma (QGP), a state in which quarks and gluons move much beyond the
size of hadrons. Heavy-ion collisions at relativistic energies are performed to create and
quantify the properties of QGP~\cite{Busza:2018rrf,Shuryak:2017aol}. Charmonia and bottomonia,
which are bound states of charm-anticharm (c$\bar{\rm c}$) or bottom-antibottom (b$\bar{\rm b}$) quarks,
respectively, are among the most sensitive probes of the characteristics of the QGP~\cite{Brambilla:2010cs}.
These bound states of heavy quarks are formed early in the heavy ion collisions and their
yields are expected to be suppressed in the medium as compared to their yields in pp
collisions ~\cite{Matsui:1986dk,Hashimoto:1986nn}.
 There have been a large number of studies on this phenomenon both theoretically and
 experimentally~\cite{Brambilla:2010cs,Schukraft:2013wba,Andronic:2015wma}
 enriching our understanding on quarkonia as probes of QGP.
 The J/$\psi$ meson was measured at the SPS, in PbPb and In-In interactions
at the centre-of-mass energy per nucleon pair $\sqrt{s_{_{NN}}}$ = 17.2 GeV~\cite{Alessandro:2004ap,Arnaldi:2007zz},
at RHIC in AuAu interactions at
$\sqrt{s_{_{NN}}}$ = 200 GeV~\cite{Adare:2011yf,Abelev:2009qaa,Tang:2011kr}, and finally at the LHC,
in PbPb collisions at $\sqrt{s_{_{NN}}}$ = 2.76 and
5.02 TeV~\cite{Aad:2010aa,Chatrchyan:2012np,Khachatryan:2016ypw,Sirunyan:2017isk,ATLAS:2016qpn,Abelev:2013ila,Adam:2016rdg,Acharya:2017tgv}.
The suppression of $J/\psi$ observed at SPS was termed as 'anomalous' $J/\psi$
suppression. It was even considered to be a hint of QGP~\cite{Alessandro:2004ap,Arnaldi:2007zz} formation.
Early theoretical calculations predicted $J/\psi$ suppression due to colour screening
in a deconfined medium which become stronger as the QGP
temperature increases~\cite{Matsui:1986dk,Digal:2001ue} but the RHIC measurements
at $\sqrt{s_{_{NN}}}$ = 200 GeV showed almost the same amount of suppression, contrary to
expectation~\cite{Brambilla:2010cs,Adare:2011yf}. These observations hint a scenario where,
at higher collision energies, the expected larger suppression is compensated by
$J/\psi$ regeneration via recombination of two independently produced
charm quarks~\cite{Andronic:2003zv,Thews:2000rj,Du:2017qkv}.

The LHC collected first PbPb collision data at $\sqrt{s_{_{NN}}} =$ 2.76 TeV at the end of 2010.
The ATLAS was the first detector to measure
the ratio of $J/\psi$ meson with Z$^{0}$ boson hinting a centrality-dependent suppression
of the yield of $J/\psi$ mesons~\cite{Aad:2010aa}.
 The $J/\psi$ measurements at high transverse momentum ($p_T>6.5$ GeV/$c$) in PbPb collisions at $\sqrt{s_{_{NN}}}$=2.76 TeV and at
$\sqrt{s_{_{NN}}}$=5.02 TeV were carried out by the CMS experiment~\cite{Chatrchyan:2012np,Khachatryan:2016ypw,Sirunyan:2017isk}.
The nuclear modification factor $R_{AA}$ of these high $p_T$ prompt $J/\psi$ decreases
with increasing centrality. The $R_{AA}$ shows a slow increase as a function of
$p_T$ after 15 GeV$/c$ and then saturates between a value 0.4 and 0.5 showing that the $J/\psi$
remains suppressed, even at very high $p_T$, up to $\sim$ 50 GeV/$c$. The ATLAS experiment
also measured $R_{AA}$ of $J/\psi$ at $\sqrt{s_{_{NN}}}$=5.02 TeV for $J/\psi$ meson having
transverse momentum above 9.0 GeV$/c$~\cite{ATLAS:2016qpn}. The amount
of suppression of $J/\psi$ mesons observed by ATLAS is similar to the suppression
observed by CMS experiment. By comparing these measurements with the STAR
results~\cite{Tang:2011kr} at RHIC, it follows that at high $p_{T}$ the suppression of
$J/\psi$ increases with collision energy.

  The ALICE experiment measured the nuclear modification factor of $J/\psi$ mesons
in the forward rapidity (2.5$<y<$4.0) range at $\sqrt{s_{_{NN}}}$=2.76 and $\sqrt{s_{_{NN}}}$=5.02
TeV~\cite{Abelev:2013ila,Adam:2016rdg} with $J/\psi$ transverse momentum starting from
0.3 GeV/$c$.
The ALICE results show that nuclear modification factor of J/$\psi$ at
low $p_{T}$ ($p_{T}<$ 12 GeV/$c$) has almost no collision centrality dependence
except in the most peripheral region where it reaches almost unity.
The $R_{AA}$ of $J/\psi$ mesons decreases substantially as a function of $p_{T}$
in the ALICE experiment~\cite{Abelev:2013ila,Adam:2016rdg}. The ALICE measurements also
give a hint for an increase of $R_{AA}$ between $\sqrt{s_{_{NN}}}$ = 2.76 and 5.02 TeV in the intermediate
$p_{T}$ region, 2 $< p_{T} < $ 6 GeV/$c$. A comparison of
the ALICE and PHENIX measurements reveal that at LHC, $J/\psi$ mesons with low $p_T$ are less
suppressed as compared to RHIC. In general, the results at LHC experiments when
compared to the RHIC measurements indicate that the data on $J/\psi$ production support a
picture where both suppression and regeneration take place in the QGP,
the two mechanisms being dominant at high and low $p_{T}$, respectively~\cite{P.ShuklaforCMS:2014vna,Kumar:2014kfa}.

In addition to $J/\psi$ mesons, the bottomonia states ($\Upsilon$(nS)) are also measured at
the LHC with very good statistical
precision~\cite{Chatrchyan:2011pe,Chatrchyan:2012lxa,Abelev:2014nua,Khachatryan:2016xxp}.
The CMS measurements at $\sqrt{s_{_{NN}}} =$2.76 TeV~\cite{Chatrchyan:2011pe,Chatrchyan:2012lxa} reveal
a clear proof of sequential suppression :  $\Upsilon$(2S) and $\Upsilon$(3S) are
more suppressed relative to the ground state $\Upsilon$(1S).   The individual $\Upsilon$ states are also found to be suppressed in
the PbPb collisions relative to the production in the pp collisions. The $\Upsilon$ nuclear
modification factor, $R_{AA}$, shows a strong dependence on collision centrality but has
weak dependence on $\Upsilon$ meson $p_{T}$ and rapidity~\cite{Khachatryan:2016xxp}.
The forward rapidity ($2.5 \leq y^{\Upsilon} \leq 4.0$) measurement of the $\Upsilon$ suppression at
ALICE~\cite{Abelev:2014nua} is found to be consistent with the midrapidity ($|y^{\Upsilon}|\,\leq 2.4$)
measurement of the $\Upsilon$ suppression at the CMS.
The CMS and ALICE collaborations have carried out the $R_{AA}$ measurement of $\Upsilon$
at $\sqrt{s_{_{NN}}} =$ 5.02 TeV with the Run II LHC PbPb
collisions~\cite{Sirunyan:2017lzi,Sirunyan:2018nsz,ALICE:Y5TeV}.
The CMS experiment measured slightly more amount of $\Upsilon$ suppression at
$\sqrt{s_{_{NN}}} =$ 5.02 TeV~\cite{Sirunyan:2017lzi,Sirunyan:2018nsz} than the suppression at
$\sqrt{s_{_{NN}}} =$ 2.76 TeV~\cite{Khachatryan:2016xxp} while the ALICE experiment observed less
suppression at $\sqrt{s_{_{NN}}} =$ 5.02 TeV than that at $\sqrt{s_{_{NN}}} =$ 2.76 TeV
in the most central PbPb collisions~\cite{Abelev:2014nua,ALICE:Y5TeV}.

 The field of quarkonia has attracted a large amount of theoretical activities
since the first prediction of $J/\psi$ suppression in heavy
ion collisions~\cite{Matsui:1986dk,Hashimoto:1986nn}.
The color screening
of q$\overline{\rm q}$ potential inside the QGP~\cite{Karsch:1987pv,Abdulsalam:2012bw}
was the first approach to study the suppression of quarkonia
in heavy ion collisions.
 In a complementary way to this static approach, $J/\psi$ suppression can also be understood
as a result of dynamical interactions with the surrounding gluons~\cite{Bhanot:1979vb,Chen:2017jje,Kumar:2014kfa}.
One of the ways to calculate the regeneration is to use the principle of detailed balance~\cite{Thews:2000rj}.
There are other effects also, namely shadowing and comover interaction~\cite{Vogt:2015uba,Ferreiro:2014bia}.
Shadowing arises as the parton distribution functions are modified inside the nucleus.
A comprehensive framework to explain the experimental data from LHC considering shadowing and
comover~\cite{Du:2017qkv,Rapp:2017chc} and viscous hydrodynamics~\cite{Krouppa:2015yoa,Krouppa:2017jlg}
has recently been formulated.

 Some of us have studied the modification of quarkonia yields due to various
processes in AA collisions in a previous work~\cite{Kumar:2014kfa}.
In that work, the gluon dissociation cross section has
been adopted from the calculations of Bhanot and Peskin~\cite{Bhanot:1979vb}
where the gluon dissociation rate has been estimated from operator product expansion
in the Coulomb approximation. Recently, Chen and He have revisited the gluon dissociation
using QCD multipole expansion for various quarkonia in QGP~\cite{Chen:2017jje}.
They reproduced the result of Bhanot and Peskin as the colour-electric
dipole (E1) transition. They have also calculated the colour magnetic
dipole (M1) transition and found its contribution to be significant
at low energies.
 In this paper, we calculate $J/\psi$ and $\Upsilon$ evolution
in a kinetic model which includes dissociation by thermal gluons
(both E1 and M1 transitions), modification
of their yields due to shadowing and due to collisions with comovers.
Regeneration by near thermal heavy quark pairs is also considered in the calculations.
We obtain the nuclear modification factor of quarkonia as a function of the
transverse momentum and centrality of the collision and compare it to the latest
experimental data from LHC at $\sqrt{s_{_{NN}}} =$ 5.02 TeV.

\section{Theoretical formulation}

The theoretical formulation towards the formation and dissociation of quarkonia have already been discussed
in Ref.~\cite{Kumar:2014kfa}. We use the same formulation here and hence discuss only some important points.

\subsection{Formation and Dissociation}

In the kinetic approach \cite{Thews:2000rj}, the evolution of the quarkonia
population $N_{Q}$ with the proper time, $\tau$, is given by the kinetic equation

\begin{equation}\label{eqkin}
{dN_{Q} \over d\tau}  =  - \lambda_D  \rho_g N_{Q} + \lambda_F {N_{q \bar{q}}^{2} \over V(\tau)}
\end{equation}
In the above equation $V(\tau)$ is the spatial volume of the QGP  and $N_{q \bar{q}}$ is
the number of initial heavy quark pairs produced in a event as a function of the centrality
defined by the number of participants $N_{\rm part}$. The $\lambda_{D}$ is the dissociation
rate and the $\lambda_{F}$ is the formation rate. Here $\rho_g$ is the density of thermal gluons.

\begin{figure}
\includegraphics[width=0.49\textwidth]{Fig1a_JPsi_SigmaDq0.pdf}
\includegraphics[width=0.49\textwidth]{Fig1b_Y1S_SigmaDq0.pdf}
\caption{(Color online) Gluon dissociation cross section of J/$\psi$ (a) and $\Upsilon$(1S) (b) as a function of gluon energy
    ($q^{0}$) in quarkonia rest frame. The total dissociation cross-section is sum of both color-electric and chromo-magnetic dipole
    transitions.}
\label{fig:SigmaDQ0}
\end{figure}

The gluon dissociation cross section of quarkonia was calculated by
Bhanot and Peskin using operator product expansion method in the color dipole
approximation~\cite{Bhanot:1979vb}. Recently, Chen and He have revisited the
gluon dissociation using QCD multipole expansion method~\cite{Chen:2017jje}.
They reproduced the result of Ref.~\cite{Bhanot:1979vb} as the color-electric dipole (E1)
transition and also calculated the color magnetic dipole (CMD-M1) transition.
The E1 transition cross-section of gluon dissociation as a function of gluon energy, $q^0$,
in the quarkonium rest frame is~\cite{Bhanot:1979vb}
\begin{equation}
\sigma^{E1}_{D}(q^{0}) = {8\pi \over 3} \, {16^2 \over 3^2} {a_0 \over m_q}  \frac{(q^0/\epsilon_0 - 1)^{3/2}} {(q^0/\epsilon_0)^5}
\end{equation}
 where $\epsilon_0$ is the quarkonia binding energy and $m_q$ is the charm/bottom quark mass
and $a_0=1/\sqrt{m_q\epsilon_0}$. Using the same procedure, the M1 transition cross-section of gluon
dissociation can be calculated as~\cite{Chen:2017jje}
\begin{equation}
\sigma^{M1}_{D}(q^{0}) = {8\pi \over 3} \, {16 \over 3 } {a_0 \over m_q^{2}}  \frac{\epsilon_0(q^0/\epsilon_0 - 1)^{3/2}} {(q^0/\epsilon_0)^3}
\end{equation}
The total dissociation cross section is given as
\begin{equation}
\sigma_D(q^0) = \sigma^{E1}_D(q^0) + \sigma^{M1}_D(q^0).
\end{equation}
The values of $\epsilon_0$ are taken as 0.64 GeV and 1.10 GeV for the ground states, $J/\psi$ and $\Upsilon$(1S),
respectively~\cite{Karsch:1987pv}. For the excited states of bottomonia, the dissociation cross sections are
used from Ref.~\cite{Chen:2017jje,Oh:2001rm}.

Figure~\ref{fig:SigmaDQ0} shows the gluon dissociation cross sections
of $J/\psi$ and $\Upsilon$(1S) as a function of gluon energy
for both color-electric and color-magnetic dipole transitions.
The color-magnetic dipole transition cross-section has similar shape as that of
electric cross-section and gives significant contribution
in low and intermediate gluon energies. Total dissociation cross section increases
with gluon energy and reaches a maximum around 0.9 GeV for
J/$\psi$ and around 1.5 GeV for $\Upsilon$(1S). At higher gluon energies, the
interaction probability decreases.
 The dissociation rate is calculated as a function of quarkonium momentum
by integrating the dissociation cross section over thermal gluon momentum
distribution $f_{g}(p_g)$ as
\begin{eqnarray}
\lambda_{D} \rho_{g}  & = & \langle \sigma v_{\rm rel} \rangle \,\rho_{g}  = \frac{g_g}{(2\pi)^{3}} \int d^{3}p_g \, f_{g}(p_g)  \, \sigma_{D}(s) v_{\rm rel}(s)  \nonumber \\
                   & = & \frac{g_g}{(2\pi)^{3}} \int dp_g 2\pi p_g^{2} f_{g}(p_g) \int d\,{\rm cos\theta}\,\sigma_{D}(s)\,v_{\rm rel}(s),
\label{Eqn:Diss}
\end{eqnarray}
where $\sigma_{D}(s) = \sigma_{D}(q^0(s))$ in terms of the square of the center of
mass energy $s$ of the quarkonium-gluon system given by
$s=M_{Q}^{2} + 2  p_g \, \sqrt{M_{Q}^2 + p^2} - 2  p_g \, p \, {\rm cos\theta}$.
Here $M_{Q}$ is the mass and $p$ is the momentum of quarkonium and $\theta$ is the angle
between the quarkonium and the gluon.
  The variables $q^0$ and $s$ are related by $q^{0} = (s-M_{Q}^{2})/(2\,M_{Q})$.
$v_{\rm rel}$ is the relative velocity between the quarkonium and the gluon~\cite{Kumar:2014kfa}.
  The formula in Eq.~\ref{Eqn:Diss} is assumed
for the most central collisions. We multiplied by a system size dependent factor
($\sqrt{N_{\rm part}/2A}$) to get the dissociation rate for other centralities.
The $J/\psi$ gluon dissociation rates as a function of temperature $T$ are shown in
Fig.~\ref{fig:DRateVsTempAndPt}(a) and as a function of $p_T$ in Fig.~\ref{fig:DRateVsTempAndPt}(b).
The dissociation rate increases with temperature as the gluon density increases.
Also, it is maximum when the quarkonium is at rest and then decreases with $p_T$.

\begin{figure}
\includegraphics[width=0.49\textwidth]{Fig2a_DRateVsT.pdf}
\includegraphics[width=0.49\textwidth]{Fig2b_DRateVsPt.pdf}
\caption{(Color online) Gluon dissociation rate of $J/\psi$ as a function of (a) temperature and
(b) $J/\psi$ transverse momentum.}
\label{fig:DRateVsTempAndPt}
\end{figure}

We can calculate the formation cross section from the dissociation cross section
using principle of detailed balance~\cite{Thews:2000rj,Thews:2005vj} as follows
\begin{equation}
\sigma_{F} = \frac{48}{36}\,\sigma_{D}(q^0)\frac{(s-M_{Q}^2)^{2}}{s(s-4m_q^{2})}.
\end{equation}
The formation rate of quarkonium at momentum {\bf p} can be written as
\begin{equation}
\frac{d\lambda_{F}}{d{\rm\bf p}} = \int \,d^{3}p_1 \,d^{3}p_2 \,\sigma_{F}(s)\, v_{\rm rel}(s)\,f_{q}(p_1)\, f_{\bar{q}} (p_2)\,\delta({\rm\bf p}-( {\rm\bf p_1} + {\rm\bf p_2} )).
\end{equation}
The variable $s$ is the square of center of mass energy between the two heavy quarks with
energy-momenta ($E_1$, ${\rm\bf p_1}$) and ($E_2$, ${\rm\bf p_2}$) with $v_{\rm rel}$ as their
relative velocity.

The functions $f_{q/\bar{q}}(p)$ are taken as normalized near-thermal distribution functions
of $q/\bar{q}$. These distributions can be described by the Tsallis
function as follows
\begin{equation}
f_{q} (p,T) = A_{n}\,\left( 1+\frac{ \sqrt{p^2+m_q^2} }{n \,T} \right)^{-n}.
\end{equation}
Here $A_n$ is the normalization factor and $n=12$ is obtained by fitting the
transverse momentum spectra of D mesons measured by CMS experiment~\cite{Sirunyan:2017xss}.
Figure~\ref{fig:Figure3_Tsallis} shows the
transverse momentum spectra of D mesons in pp and PbPb collisions
at $\sqrt{s_{\rm NN}}=5.02$ TeV measured by CMS experiment. The spectra are fitted
by Tsallis function with $n=6.9$ for pp and $n=12$ for PbPb collisions.
\begin{figure}
\begin{center}
  \includegraphics[width=0.60\textwidth]{Fig3_DMeson_502TeV.pdf}
  \caption{The transverse momentum spectra of D mesons in pp and PbPb collisions
    at $\sqrt{s_{\rm NN}}=5.02$ TeV measured by CMS experiment~\cite{Sirunyan:2017xss}.
    The spectra are fitted by Tsallis function with $n=6.9$ for pp and $n=12$ for PbPb collisions.
}
\label{fig:Figure3_Tsallis}
\end{center}
\end{figure}

Figure~\ref{fig:ForRateVsTempAndPt}(a) shows the behaviour of the formation rate as a function
of temperature at different values of $J/\psi$ meson $p_{T}$.
Figure~\ref{fig:ForRateVsTempAndPt}(b) shows the same as a function
of $J/\psi$ meson $p_T$ at different temperatures.
It shows that the $J/\psi$ generated from uncorrelated heavy quark pairs has
softer $p_{T}$ distribution than that of $J/\psi$'s coming from the initial hard scatterings.
Thus the effect of recombination will be important at low and intermediate $p_T$.

\begin{figure}
\includegraphics[width=0.49\textwidth]{Fig4a_FormRateVsT.pdf}
\includegraphics[width=0.49\textwidth]{Fig4b_FormRateVsPt.pdf}
\caption{(Color online) Formation rate of  $J/\psi$ as a function of (a) temperature and
(b) $J/\psi$ transverse momentum.}
\label{fig:ForRateVsTempAndPt}
\end{figure}

\subsection{Nuclear modification}

The nuclear modification factor ($R_{AA}$) as a function of $p_T$ can be obtained as

\begin{equation}
  R_{AA}(p_T)= \frac{ \Sigma_{\rm Centrality} \, N_Q(p_T,N_{\rm part})}
                    { \Sigma_{\rm Centrality} \, N_{\rm coll} \, N_{Q}^{pp}(p_T)}.
\label{raa}
\end{equation}
The sum over the events is performed over the measured centrality range in the
experiment and N$_{\rm coll}$ is the number of binary collisions for the centrality bins
used in experiments.
 The $R_{AA}$ as a function of collision centrality is obtained as

\begin{equation}
  R_{AA}(N_{\rm part}) = \frac{\int_{p_{T\,\rm cut}} N_{Q}(p_T, N_{\rm part}) \, p_T dp_T}
                  {\int_{p_{T\,\rm cut}} N_{\rm coll} \,\, N^{\rm pp}_{Q}(p_T) \, p_T dp_T}
\label{raa2}
\end{equation}
Here $p_{T~{\rm cut}}$ defines the $p_T$ range for acceptance of the experiment.
{\color{black}The shape of quarkonia $p_{T}$ distribution in pp collisions is obtained from
  PYTHIA 6.424~\cite{Sjostrand:2006za}. The PYTHIA generator gives good description of quarkonia data
  at LHC energy~\cite{Khachatryan:2010yr}. The number of Q$\bar{Q}$ pairs in different centrality classes
  are obtained by the N$_{\rm coll}$ scaling. The N$_{\rm coll}$ values in different centrality classes
  are calculated using the Glauber model~\cite{Loizides:2017ack}.
}

The evolution of the system for each centrality range is governed by an isentropic
cylindrical expansion ($s(T)\,V(\tau)= s(T_0)\,V(\tau_0)$) with prescription given
in Ref.~\cite{Kumar:2014kfa}.
The equation of state (EOS) obtained by Lattice QCD and by hadronic resonances is
used~\cite{Huovinen:2009yb}. The radius $R$ for a given centrality
with number of participant $N_{\rm part}$ is obtained as $R(N_{\rm part}) = R_{A} \sqrt{N_{\rm part}/2A}$,
where $R_{A}$ is radius of the nucleus.
The initial entropy density, $s(\tau_0)|_{0-5\%}$, for 0-5\% centrality is
\begin{eqnarray}
s(\tau_0)|_{0-5\%}  = {a_{m} \over V(\tau_0)|_{0-5\%}}   \left(\frac{dN}{d\eta}\right)_{0-5\%} .
\label{TempVsMult}
\end{eqnarray}
Here $a_m$ (= 5) is a constant connecting the total entropy to the final hadron
multiplicity $dN/d\eta$ obtained from hydrodynamic calculations~\cite{Shuryak:1992wc}.
The initial temperature, $T_0$, in the 0-5$\%$ most central collisions is estimated
from the total multiplicity in the given rapidity, assuming that the initial time is
$\tau_0 = 0.3$ fm/$c$. The total multiplicity in a given rapidity window is
1.5 times the measured charged particle multiplicity in PbPb collisions at 5.02 TeV.
With the lattice EOS, at midrapidity, with $(dN_{\rm ch}/d\eta)_{0-5\%} = 1943$~\cite{Adam:2015ptt},
we find $T_0 = 0.516$ GeV. Likewise, at forward rapidity~\cite{Adam:2016ddh}; $T_0$ is
0.487 GeV. The freeze out temperature is taken to be $T_f=0.140$ GeV.

\subsection{Hadronic comovers}

The suppression of quarkonia caused due to comoving pions is obtained by folding the quarkonium-pion
dissociation cross section $\sigma_{\pi Q}$ over thermal pion distributions \cite{Vogt:1988fj}.
The  cross section of pion-quarkonia is calculated by convoluting the gluon-quarkonia cross section $\sigma_D$
over the gluon distribution obtained inside the pion~\cite{Arleo:2001mp},
\begin{equation}
\sigma_{\pi Q} (p_{\pi}) = {p_+^2 \over 2(p_\pi^2 - m_\pi^2)} \int_0^1 \, dx \, G(x) \, \sigma_D(xp_+/\sqrt {2}),
\end{equation}
where $p_+ = (p_\pi + \sqrt{p_\pi^2-m_\pi^2})/\sqrt{2}$. The term  $G(x)$ is the gluon distribution inside a pion which
can be given by the GRV parameterization~\cite{Gluck:1991ey}.
The comover cross section  is expected  to be small  at LHC energies~\cite{Lourenco:2008sk}.

\subsection{Experimental data to fit}

\begin{table*}[b]
\caption{Total $\rm{c\overline{c}}$ production cross-section measured by ALICE experiment in pp collisions at LHC.}
\begin{tabular}{l|l|l|l}
\hline
\hline
  $\sqrt{s}$(TeV)           &$\sigma^{c\bar{c}}\pm {\rm stat.}\pm {\rm syst.}$(mb)            &Experiment      &Ref.  \\
\hline
 2.76 TeV                   &$4.8\pm0.8^{+1.0}_{-1.3}$           & ALICE     &\cite{Abelev:2012vra}               \\
 7 TeV                      &$8.5\pm0.5^{+1.0}_{-2.4}$           & ALICE     &\cite{Abelev:2012vra}               \\
 \hline
\hline
\end{tabular}
\label{CCBarCrossExp}
\end{table*}

\begin{table}[t]
\caption[]{Total $\rm{b\overline{b}}$ production cross-section measured by ALICE experiment in pp collisions at LHC.}
\label{BBBarCrossExp}
\begin{tabular}{l|l|l|l}
\hline
\hline
  $\sqrt{s}$(TeV)           &$\sigma^{\rm{b\overline{b}}}\pm {\rm stat.}\pm {\rm syst.}$($\mu$b)            &Experiment      &Ref.  \\
\hline
 2.76 TeV                   &$130\pm15.1^{+42.1}_{-49.8}$           & ALICE     &\cite{Abelev:2014hla}               \\
 7 TeV                      &$282\pm74^{+58}_{-68}$               & ALICE     &\cite{Abelev:2012gx}               \\
\hline
\hline
\end{tabular}
\end{table}

The total charm and total bottom production cross-sections are measured by different experiments at
LHC~\cite{Abelev:2012vra,Abelev:2014hla,Abelev:2012gx}.
Table~\ref{CCBarCrossExp} and Table~\ref{BBBarCrossExp} show the values of total $\rm{c\overline{c}}$
and $\rm{b\overline{b}}$ production cross-section measured in pp collisions at LHC.
We use these values to estimate the
total heavy quark production cross section at $\sqrt{s_{_{NN}}} =$ 5.02 TeV.
The quarkonium production cross sections are calculated from the measured heavy
quark production cross-section using the energy independent factors (0.00526 for J/$\psi$ and 0.002 for $\Upsilon$)
obtained from the color evaporation model~\cite{Kumar:2014kfa,Nelson:2012bc,Vogt:2012vr}.
The cold nuclear matter (CNM) effects i.e. the modifications of the parton distribution
functions (nPDF) in PbPb collisions is calculated using the central EPS09 NLO
parameter set~\cite{Eskola:2009uj}.
The uncertainty in cold matter effects is obtained by adding the EPS09 NLO uncertainties in quadrature.

The production cross sections for heavy flavor and quarkonia at $\sqrt{s_{_{NN}}} =$ 5.02
TeV are given in Table~\ref{Tab:NLOcross}.  The yields in a minimum bias
PbPb event is obtained from the per nucleon cross
section, $\sigma_{\rm PbPb}$ as
\begin{eqnarray}
  N = \frac{A^2 \sigma_{\rm PbPb}}{\sigma_{\rm PbPb}^{\rm tot}} \, \, .
\end{eqnarray}
At 5.02 TeV, the total PbPb cross section, $\sigma_{\rm PbPb}^{\rm tot}$,
is 7.7 b~\cite{Loizides:2017ack}.

\begin{table}
\caption[]{ Heavy quark and quarkonia production  cross sections per nucleon pair at
  $\sqrt{s_{_{_{NN}}}}= 5.02$ TeV. The quantity $N^{\rm PbPb}$ gives the initial number of
  heavy quark pair/quarkonia in one PbPb event.}
\label{Tab:NLOcross}
\begin{tabular}{l|l|l|l|l}
\hline
\hline
             &$ c \overline c$            &J/$\psi$                      & $ b \overline b$                    & $\Upsilon$   \\
\hline
$\sigma_{pp}$  &$6.754^{+1.195}_{-2.015}$ mb     & $35.32^{+6.247}_{-10.537}~\mu$b     & $210.30^{+70.769}_{-77.640}~\mu$b            & $0.4206^{+0.142}_{-0.155}~\mu$b  \\

&&&&\\
$\sigma_{\rm PbPb}$ &$4.669^{+0.826}_{-1.393}$ mb    &$24.56^{+4.344}_{-7.327}~\mu$b    & $179.30^{+60.337}_{-66.195}~\mu$b             &$0.3586^{+0.121}_{-0.132}~\mu$b  \\
&&&&\\

$N^{\rm PbPb}$     &$26.23^{+4.639}_{-7.825}$       & $0.1381^{+0.024}_{-0.041}$         & $1.007^{+0.339}_{-0.372}$                  & $0.0020^{+0.0007}_{-0.0007}$       \\

\hline
\hline
\end{tabular}
\end{table}

\section{Results and discussions}

\begin{figure}
\begin{minipage}{1.0\linewidth}
\centering
{\includegraphics[width=0.49\textwidth]{Fig5a_CMS_RAAPt_Shade.pdf}}
{\includegraphics[width=0.49\textwidth]{Fig5b_ATLAS_RAAPt_Shade.pdf}}
\end{minipage}%
\ \\
\centering
\begin{minipage}{0.5\linewidth}
\centering
{\includegraphics[width=1.0\textwidth]{Fig5c_ALICE_RAAPt_Shade.pdf}}
\end{minipage}%
\caption{(Color online)
Calculated nuclear modification factor ($R_{AA}$) as a function of $J/\psi$
transverse momentum compared with (a) CMS, (b) ATLAS and (c) ALICE
measurements~\cite{Sirunyan:2017isk,ATLAS:2016qpn,Adam:2016rdg}.
The global uncertainty in $R_{AA}$ is shown as a band around the line at 1.
}
\label{fig:JPsiRaaVsPt}
\end{figure}

\begin{figure}
\begin{minipage}{1.0\linewidth}
\centering
{\includegraphics[width=0.49\textwidth]{Fig6a_CMS_RAANPart_Shade.pdf}}
{\includegraphics[width=0.49\textwidth]{Fig6b_ATLAS_RAANPart_Shade.pdf}}
\end{minipage}%
\ \\
\centering
\begin{minipage}{0.5\linewidth}
\centering
{\includegraphics[width=1.0\textwidth]{Fig6c_ALICE_RAANPart_Shade.pdf}}
\end{minipage}%
\caption{(Color online) Calculated nuclear modification factor ($R_{AA}$) of $J/\psi$ as a function
  of centrality of collisions, compared with (a) CMS, (b) ATLAS and (c) ALICE
  measurements~\cite{Sirunyan:2017isk,ATLAS:2016qpn,Adam:2016rdg}.
The global uncertainty in $R_{AA}$ is shown as a band around the line at 1.
}
\label{fig:JPsiRaaVsNPart}
\end{figure}

\begin{figure}
\includegraphics[width=0.49\textwidth]{Fig7a_Y1S_CMS_RAAPt_Shade.pdf}
\includegraphics[width=0.49\textwidth]{Fig7b_Y2S_CMS_RAAPt_Shade.pdf}
\caption{(Color online) Calculated nuclear modification factor ($R_{AA}$) of (a) $\Upsilon$(1S) and
  (b) $\Upsilon$(2S) as a function of $p_{T}$
  compared with CMS measurements~\cite{Sirunyan:2018nsz}.
The global uncertainty in $R_{AA}$ is shown as a band around the line at 1.
}
\label{fig:UpsilonRaaPtCMS}
\end{figure}

\begin{figure}
\includegraphics[width=0.49\textwidth]{Fig8a_ALICE_Y1SRAAPt_Shade.pdf}
\includegraphics[width=0.49\textwidth]{Fig8b_ALICE_Y2SRAAPt_Shade.pdf}
\caption{(Color online) Calculated nuclear modification factor ($R_{AA}$) of (a) $\Upsilon$(1S) and
  (b) $\Upsilon$(2S) as a function of $p_{T}$ in the kinematic range of ALICE detector at LHC ~\cite{ALICE:Y5TeV}.
  The global uncertainty in $R_{AA}$ is shown as a band around the line at 1.
}

\label{fig:UpsilonRaaPtALICE}
\end{figure}

\begin{figure}
\includegraphics[width=0.49\textwidth]{Fig9a_CMS_Y1SRAANPart_Shade.pdf}
\includegraphics[width=0.49\textwidth]{Fig9b_CMS_Y2SRAANPart_Shade.pdf}
\caption{(Color online) Calculated nuclear modification factor ($R_{AA}$) of
  (a) $\Upsilon$(1S) and (b) $\Upsilon$(2S) as a function of centrality of the
  collisions compared with the CMS measurements~\cite{Sirunyan:2018nsz}
  The global uncertainty in $R_{AA}$ is shown as a band around the line at 1.
}
\label{fig:UpsilonRaaNPartCMS}
\end{figure}

Figure~\ref{fig:JPsiRaaVsPt}(a) and (b) show the estimations of different contributions
to the nuclear modification factor, $R_{AA}$, for the J/$\psi$ meson as a function of $p_T$
along with the mid rapidity and high $p_T$ measurements from
CMS~\cite{Sirunyan:2017isk} and ATLAS~\cite{ATLAS:2016qpn} experiments respectively.
Figure~\ref{fig:JPsiRaaVsPt}(c) shows the same for the low $p_T$ and forward
rapidity compared with the measurement by ALICE experiment~\cite{Adam:2016rdg}.
At low $p_T$, regeneration of J/$\psi$ gives dominant
contribution which overcomes the strong suppression by gluon dissociation.
This looks to be the reason for the increase of the $R_{AA}$ of J/$\psi$  around $p_T\approx$ 2 GeV/c.
The suppression due to gluon dissociation is substantial at low $p_T$ and reduces as
we move to higher $p_T$. At low and intermediate $p_T$, both regeneration and dissociation, due
to the presence of QGP are effective.  The high $p_T$
suppression ($p_T > 10$  GeV/$c$) of $J/\psi$ measured by CMS is greater than the
suppression caused by gluon dissociation in the QGP. The suppression measured
by CMS and ATLAS is at high $p_T$ values greater than the heavy quark mass,
and thus here the energy loss from initial partonic scatterings might play a
crucial role as it does for open heavy flavour.
{\color{black}
The feed-down contributions for the J/$\psi$ are not very
large. At the LHC energies, around 80\% J/$\psi$ are from the directly produced hard scattering~\cite{Lansberg:2019adr}.
Moreover, the gluon dissociation cross section for the excited charmonia states are not reliable. Thus we chose
not to consider the feed-down from the higher states for the J/$\psi$.
}
The band around the total $R_{AA}$ includes all the uncertainties in the dissociation and regeneration
processes as discussed in the next paragraph. The band around CNM effects is shown
separately and is not included in the total uncertainty band for a better display.
The values of gluon-quarkonia cross section ($\sigma_{D}$) and the initial temperature $T_{0}$
can have uncertainties and will affect the results.
The value of $\sigma_D$ is varied by $\pm$50\% around the calculated value to obtain
the uncertainty in $R_{AA}$.  The initial temperature is calculated using measured
charged particle density and nominal value of $\tau_0 =$ 0.3 fm/$c$. The value of $\tau_0$
is varied  in the range 0.1 $<\,\tau_0\,<$ 0.6 fm/$c$ which corresponds to the
variation in the initial temperature from +45 \% to - 20 \%.
The uncertainty in the charm pair cross section is also considered as a source while
obtaining the contribution due to regeneration.
The total uncertainty is obtained by adding in quadrature all individual uncertainties.

The calculation of $R_{AA}$ of J/$\psi$ is also made as a function of collision
centrality (system size).
Figure~\ref{fig:JPsiRaaVsNPart} shows calculations of different contributions to the J/$\psi$
nuclear modification factor as a function of system size, along with the measurements
from CMS in (a), ATLAS in (b) and ALICE in (c)~\cite{Sirunyan:2017isk,ATLAS:2016qpn,Adam:2016rdg}.
The figure shows that the suppression of J/$\psi$ due to
QGP increases when the system size grows. The contribution from regeneration process
is minimum in the high $p_T$ range ($p_T~>9$ GeV/c) for ATLAS case
as shown in Fig.~\ref{fig:JPsiRaaVsNPart}(b) and maximum in the low $p_T$ range
for the ALICE case shown in Fig.~\ref{fig:JPsiRaaVsNPart}(c). For low $p_T$ ALICE measurement,
the nuclear modification factor is almost flat from mid to central collisions because
the regeneration compensates the gluon dissociation, a trend which is well-reproduced
by our calculations.
The CMS centrality dependence of the $R_{AA}$ of $J/\psi$ given in
Fig.~\ref{fig:JPsiRaaVsNPart}(a) is reasonably well described by the model since
the contribution of the CMS data comes from $J/\psi$ with 6.5 $<\,p_T<\,$10 GeV/$c$
where the suppression due to gluon dissociation dominates.
For the case of ATLAS data, the model reproduces the shape of centrality dependence
of the $R_{AA}$ of $J/\psi$ observed in the data.

\begin{figure}
\includegraphics[width=0.49\textwidth]{Fig10a_Y1S_ALICE_RAANPart_Shade.pdf}
\includegraphics[width=0.49\textwidth]{Fig10b_Y2S_ALICE_RAANPart_Shade.pdf}
\caption{(Color online) Calculated nuclear modification factor ($R_{AA}$) of
  (a) $\Upsilon$(1S) and (b) $\Upsilon$(2S) as a function of centrality of
  the collisions compared with the ALICE measurement~\cite{ALICE:Y5TeV}.
The global uncertainty in $R_{AA}$ is shown as a band around the line at 1.
}
\label{fig:UpsilonRaaNPartALICE}
\end{figure}

Figure~\ref{fig:UpsilonRaaPtCMS}(a) and (b) show the calculations of contributions to
the nuclear modification factor, $R_{AA}$, for the $\Upsilon$(1S) and $\Upsilon$(2S)
respectively as a function of $p_T$ compared with the mid rapidity measurements from
CMS~\cite{Sirunyan:2018nsz}.
The gluon dissociation mechanism combined with the pion dissociation and shadowing
corrections gives good description of data in mid $p_{T}$ range ($p_{T}\approx$ 5-10 GeV/c)
for both $\Upsilon$(1S) and $\Upsilon$(2S).
The contribution from the regenerated $\Upsilon$s is negligible even at LHC energies.
Our calculations under-predict the suppression observed at the highest measured
$p_{T}$ for $\Upsilon$(1S) and $\Upsilon$(2S) which is similar for the case
of J/$\psi$.
The states $\Upsilon$(1S) and $\Upsilon$(2S) also have
feed-down contributions from decays of higher b$\bar{\rm b}$ bound states.
The nuclear modification factor, $R_{AA}$ is obtained taking into account the feed-down corrections as follows
  \begin{equation}
    R_{AA}^{\Upsilon(3S)} = R_{AA}^{\Upsilon(3S)}\\
  \end{equation}
  \begin{equation}
    R_{AA}^{\Upsilon(2S)} = f_1 R_{AA}^{\Upsilon(2S)} +  f_2 R_{AA}^{\Upsilon(3S)} \\
  \end{equation}
   \begin{equation}
    R_{AA}^{\Upsilon(1S)} = g_1 R_{AA}^{\Upsilon(1S)} +  g_2 R_{AA}^{\chi_b(1P)} + g_3 R_{AA}^{\Upsilon(2S)} + g_4 R_{AA}^{\Upsilon(3S)}\\
  \end{equation}
The factors $f$’s and $g$’s are obtained from CDF measurement~\cite{Affolder:1999wm}.
The values of $g_1$, $g_2$, $g_3$ and $g_4$ are 0.509, 0.27, 0.107
and 0.113 respectively. Here $g_4$ is assumed to be the combined fraction of
$\Upsilon$(3S) and $\chi_b$(2P).
The values of $f_1$ and $f_2$ are taken as 0.50~\cite{Strickland:2011aa}.

Figure~\ref{fig:UpsilonRaaPtALICE}(a) and (b) show the model
prediction of the nuclear modification factor, $R_{AA}$, for the $\Upsilon$(1S)
and $\Upsilon$(2S) respectively as a function of $p_T$ in the kinematic range
covered by ALICE detector. The ALICE data~\cite{ALICE:Y5TeV} is well described by our model.

Figure~\ref{fig:UpsilonRaaNPartCMS}(a) depicts the calculated
centrality dependence of the $\Upsilon$(1S) nuclear
modification factor, along with the midrapidity data from CMS~\cite{Sirunyan:2018nsz}.
Our calculations combined with the pion dissociation and shadowing corrections
gives very good description of the measured data. Figure~\ref{fig:UpsilonRaaNPartCMS}(b)
shows the same for the $\Upsilon$(2S) along with the midrapidity
CMS measurements. The suppression of the excited $\Upsilon$(2S) states
is also well described by our model. As stated earlier, the effect of regeneration is
negligible for $\Upsilon$ states.

Figure~\ref{fig:UpsilonRaaNPartALICE}(a) shows the forward rapidity ALICE
measurement of the $\Upsilon$(1S) nuclear modification factor~\cite{ALICE:Y5TeV}
along with our calculations. The suppression due to thermal gluon dissociation
describes the measured data after including the comover and shadowing corrections.
Figure~\ref{fig:UpsilonRaaNPartALICE}(b) shows the calculations for the
$\Upsilon$(2S) nuclear modification factor in ALICE detector kinematic range.
The suppression due to thermal gluon dissociation describes the
ALICE measurements after including the comover and shadowing corrections.

\section{Summary}
 We presented detailed calculations of the $J/\psi$ and $\Upsilon$
production and the modifications of their yields in PbPb collisions at $\sqrt{s_{_{NN}}} =$ 5.02 TeV.
A kinetic model is employed which incorporates quarkonia suppression inside QGP, suppression
due to hadronic comovers and regeneration from heavy quark pairs.
The behaviour of the dissociation and formation rates are studied as a function of
transverse momentum and medium temperature.
The nuclear modification factors for both $J/\psi$ and $\Upsilon$ are obtained
as a function of system size and transverse momentum and have been compared to the measurements
in PbPb collisions at $\sqrt{s_{_{NN}}} =$ 5.02 TeV.
It is found that  regeneration of $J/\psi$ is the dominant process at low $p_T$. As a result
the $J/\psi$ production is found to be enhanced in the ALICE low $p_T$ data.
In the same $p_T$ range gluon dissociation is also substantial however it becomes small
as we move to high $p_T$.
Both  these processes affect the yields of quarkonia in a QGP medium  at low and
intermediate $p_T$.
 The high $p_{T}$ suppression ($p_T > 10$  GeV/$c$) of $J/\psi$ measured by CMS and
ATLAS is more than the suppression expected due to gluon dissociation in QGP.
The energy loss from initial partonic scatterings might play a crucial role in this region
as it does for open heavy flavour.
As the system size grows  $J/\psi$'s are increasingly suppressed.  The nuclear
modification factor at low $p_T$  (ALICE case) as a function of centrality remains
flat since the increased suppression is compensated by regenerated  $J/\psi$'s
as the system size grows.
 We could reproduce the centrality dependence of $R_{AA}$ for high $p_T$ $J/\psi$'s reasonably
 well. The $p_{T}$ and centrality dependence of suppression of $\Upsilon$ states are well reproduced
 by the model.

\section{Acknowledgement}
Authors thank Board of Research in Nuclear Sciences (BRNS) for support.
AB thanks Alexander von Humboldt (AvH) foundation and Federal Ministry
of Education and Research (Gernamy) for support through Research Group
Linkage Programme.

\section*{References}


\section{Introduction}
\label{sec:intro}

The proliferation of the Internet of Things (IoT) technology allows the concept of Smart Cities become feasible, and these IoT based applications have greatly improved the lives of residents through enhancing health, safety, and convenience. However, the first generation of IoT deployment (IoT 1.0) for smart city has come across several challenges that prevent wider adoption of the IoT in smart communities and smart cities \cite{ramachandran2019towards}. To enhance the adoption of the IoT in smart communities and smart cities, a marketplace architecture and middleware platform, called the Intelligent IoT Integrator (I3), is proposed to address key problems in IoT marketplace \cite{krishnamachari2018i3}. To build smart communities, I3 is aimed to build a ''data rivers'' that allow data streams from different entities to be merged together and analyzed, processed, and acted upon as needed to support a diverse set of applications \cite{krishnamachari2018i3}. In a commercial aspect, I3 tries to build a marketplace that allows a large volume of devices owners to sell data streams, while different applications can buy one or more data streams used for their data analytic services.

While the I3 system facilitates the easy exchange of real-time IoT data streams between devices owners and third-party applications through marketplace, it also brings new concerns on architecture and security issues. The I3 system is deployed in a distributed network environment that includes a large number of IoT devices with high heterogeneity and dynamics and different application vendors using non-standard development technology. It requires a scalable, flexible and lightweight system architecture that supports fast development and easy deployment among participants. In addition, a fundamental principle embodied in the vision of I3 is that device and data owners should have the right to decide when they will share their data (or access to their devices, in case of actuators), with whom, and at what price \cite{krishnamachari2018i3}. Thus, flows of data are mediated through end-to-end agreements between devices owners and application developers, which means that enforcing a centralized management policy on I3 becomes difficult and inefficient. Furthermore, a scalable, flexible and fine-grained data sharing and access control mechanism is needed to protect data exchange and payment among entities.

An earlier work has explored how a decentralized data marketplace (DDM) could be created using blockchain and other distributed ledger technologies \cite{ramachandran2018towards}. The possible benefits of such a decentralized architecture are highlighted and the preliminary work shows how key elements of such a decentralized marketplace could be implemented using smart contracts. Other marketplaces in this space include Ocean Protocol \cite{oceanpocotol} and IOTA Data Marketplace \cite{iota}. One of the challenges in implementing an efficient DDM lies in the management of system complexity. With many possible moving parts and components, the decentralized marketplace system could get fragmented or hard to scale. On the other hand, keeping too rigid a framework could also hurt scalability.

In this position paper, a BLockchain-ENabled Secure Microservices for DDM (BlendSM-DDM) is proposed to secure data exchange among different devices owner and application developers. Leveraging the fine-granularity and loose-coupling features of the microservices architecture \cite{butzin2016microservices, datta2018next}, the BlendSM-DDM decouples the DDM applications and security functionality into multiple containerized microservices that are computationally affordable to each individual devices. The security mechanism is implemented as separated microservices that are built on the smart contract. The distributed microservices could cooperate with each other as a service pool to perform complicated data analytical missions and security enforcement without hurting the flexibility and scalability.

The remainder of this paper is organized as follows: Section \ref{sec:relatedwork} provides the background knowledge and reviews the state of the arts research in blockchain and microservice for IoT-based system. Section \ref{sec:architecture} illustrates the details of core security features introduced by the proposed BlendSM-DDM. Section \ref{sec:conclusion} concludes this paper with our ongoing work.

\section{Related Works}
\label{sec:relatedwork}

\subsection{Microservices in IoT}
A service-oriented architecture (SOA) is widely adopted in the development of application software in an IoT environment \cite{butzin2016microservices}. The traditional SOA utilizes a monolithic architecture that constitutes different software features in a single interconnected and interdependent application and database. Owing to the tightly coupled dependence among functions and components, such a monolithic framework is difficult to adapt to new requirements in an IoT-enabled system, such as scalability, service extensibility, data privacy, and cross-platform interoperability \cite{datta2018next}. The \textit{microservices architecture} allows functional units of an application to work independently with a loose coupling though encapsulating a minimal functional software module as a microservice, which can be individually developed and deployed. finally, multiple decentralized individual microservices cooperatively perform the functions of complex systems.
The flexibility of microservices enables continuous, efficient, and independent deployment of application function units. As two most significant features of the microservices architecture, \textit{fine granularity} means each of the microservices can be developed in different frameworks and with minimal development resources, while \textit{loose coupling} implies that functions of microservices and its components are independent of each other's deployment and development \cite{yu2018survey}.

Thanks to the fine-granularity and loose-coupling properties, the microservices architecture has been investigated in many smart solutions to enhance the scalability and security of IoT-based applications. The IoT systems are advancing from ``things''-oriented ecosystem to a widely and finely distributed microservices-oriented-ecosystem \cite{datta2018next}. An Intelligent Transportation Systems (ITS) that incorporates and combines the IoT approaches using the serverless microservices architecture has been designed and implemented to help the transportation planning for the Bus Rapid Transit (BRT) systems \cite{herrera2018smart}. To enable a more scalable and decentralized solution for advanced video stream analysis for large volumes of distributed edge devices, a system design of a robust smart surveillance systems was proposed based on microservices architecture and blockchain technology \cite{nagothu2018microservice, nikouei2019decentralized, xu2019blendmas}. It aims at offering a scalable, decentralized and fine-grained access control solution for smart public safety.

\subsection{Blockchain and Smart Contract}
As a fundamental technology of Bitcoin \cite{nakamoto2008bitcoin}, \textit{blockchain} initially was used to promote a new cryptocurrency that performs commercial transactions among independent entities without relying on a centralized authority, like banks or government agencies. Essentially, the blockchain is a public ledger based on consensus rules to provide a verifiable, append-only chained data structure of transactions. Thanks to the decentralized architecture that does not rely on a centralized authority, blockchain allows the data to be stored and updated distributively. The transactions are approved by miners and recorded in the time-stamped blocks, where each block is identified by a cryptographic hash and chained to preceding blocks in a chronological order. In a blockchain network, a \textit{consensus mechanism} is enforced on a large amount of distributed nodes called miners to maintain the sanctity of the data recorded on the blocks. Thanks to the “trustless” proof mechanism running on miners across the network, users can trust the system of the public ledger stored worldwide on many different decentralized nodes maintained by ''miner-accountants'', as opposed to having to establish and maintain trust with a transaction counter-party or a third-party intermediary \cite{swan2015blockchain}. Thus, blockchain is an ideal decentralized architecture to ensure distributed transactions among all participants in a trustless environment, like edge-based IoT networks.

Emerging from the intelligent property, a \textit{smart contract} allows users to achieve agreements among parties through a blockchain network. By using cryptographic and security mechanisms, a smart contract combines protocols with user interfaces to formalize and secure relationships over computer networks \cite{szabo1997formalizing}. A smart contract includes a collection of pre-defined instructions and data that have been saved at a specific address of blockchain as a Merkle hash tree, which is a constructed bottom-to-up binary tree data structure. Through exposing public functions or application binary interfaces (ABIs), a smart contract interacts with users to offer the predefined business logic or contract agreement.

Other distributed data sharing platforms in this space includes Trinity \cite{ramachandran2019trinity}, which is a distributed pub-sub data broker with support for blockchain-based immutable data storage, Pub-Pay-Sub Protocol \cite{ramachandran2019publish}, which is a publish-subscribe broker with built-in support for micropayments, and SDPP \cite{radhakrishnan2018streaming}, which is a peer-to-peer streaming data protocol with micropayments support.

The blockchain and smart contract enabled security mechanism for applications has been a hot topic and some efforts have been reported recently, for example, smart surveillance system \cite{nagothu2018microservice, nikouei2018realtime, xu2019blendmas}, social credit system \cite{xu2018constructing}, space situation awareness \cite{xu2018exploration}, biomedical imaging data processing \cite{xu2019decentralized}, identification authentication \cite{hammi2018bubbles} and access control \cite{xu2018blendcac, xu2018smartcac}. Blockchain and smart contract together are promising to provide a solution to enable a secured data sharing and access authorization in decentralized IoT systems.

\section{System Architecture of BlendSM-DDM}
\label{sec:architecture}
In the vision of I3, the system is oriented for device owners and data users, and the data exchanged payment are mainly between sellers and buyers. To improve the economically self-sustaining business ecosystem, blockchain is a very promising solution that will enrich the I3 system by providing a trusted sharing service, where data is reliable, immutable and auditable. The virtualization technology, like virtual machines (VMs) or containers, is platform independent and could provide resource abstraction and isolation features, they are ideal for system architecture design to address the heterogeneity challenge in IoT system. Compared to VMs, containers are more lightweight and flexible with OS-level isolation, so that is an ideal selection for edge computing.

\begin{figure} [t]
\begin{center}
\begin{tabular}{c}
\includegraphics[height=11.0 cm]{fig1_Blend_I3_architecture.png}
\end{tabular}
\end{center}
\caption[example] { \label{fig:I3_architecture} Illustration of BlendSM-DDM System Architecture.}
\vspace{-10pt}
\end{figure}

Inspired by blockchain and containerized microservices technology, BlendSM-DDM is proposed as a decentralized security services architecture for DDM system, which is shown in Fig. \ref{fig:I3_architecture}. Trinity \cite{ramachandran2019trinity} distributes the pub-sub broker using the blockchain technology to prevent byzantine faults. A similar approach could be adopted to distribute the microservices among consensus nodes. A permissioned blockchain infrastructure running on a decentralized peer-to-peer network is managed by DDM administrators, and it could provide a reliable, traceable and immutable distributed ledger for trusted data sharing services in DDM system. Utilizing the containerized microservices architecture, the key elements and security features of DDM system are decoupled into multiple microservices and deployed on distributed computing hosts. Each microservice exposes REST APIs to accept service requests, and API gateway could be managed by either DDM administrator or service providers. Those decentralized microservices work as a blockchain-enabled services pool to offer a scalable, flexible and lightweight security and privacy protection mechanism for data marketplace.

\subsection{Permissioned Blockchain Network}
The permissioned blockchain network in the BlendSM-DDM requires that only registered entities are authorized access to the network, so that they can interact with the blockchain and perform blockchain-based services, such as mining block, sending transactions and deploying smart contracts. The upper part of Fig. \ref{fig:I3_architecture} illustrates the identity authentication process to enroll a new participant in the permissioned blockchain network. Given computation capacity of host devices and management policy, all participants are categorized as validator (miner) or node (non-miner). The validators are able to cooperatively executing consensus algorithm to maintain the integrity and immutability of distributed ledger. While nodes are only capable of launching transactions to synchronize data with the blockchain.

To join the permissioned blockchain network, the participants must send joining requests to an oracle for registration and identity authentication. The administrator of DDM system could work as such oracle, and maintain a global identity profile and authorization policies for the permissioned blockchain network management. After approved by the oracle, the new entity's node information will be added to a global static node record, and will be broadcasted to all certificated participants in the permissioned blockchain network accordingly. The membership revocation occurs when an entity explicitly launches a leaving request or the oracle implicitly rule out any misbehaved nodes.

Compared to a public blockchain network that suffers from performance and security issues, such as low throughput, poor privacy and lack of governance, the permissioned blockchain can enforce a more efficient but lightweight consensus mechanism to improve the transaction rate and data storage capacity. Furthermore, through limiting participants with clearly defined privacy policies, the sensitive data could be only visible to the authorized participants but not to the public. Therefore, the permissioned blockchain network also enable a more strict privacy policy than a public blockchain network does.

\subsection{Blockchain-enabled Security Services}

The blockchain-enabled security services layer, as shown in lower part of Fig. \ref{fig:I3_architecture}, acts as a fundamental microservices oriented infrastructure to support business logic functions of the DDM system as well as security mechanism. The key elements and operations are described below.

    \noindent \emph{1. Data Pub-sub Payment}: The seller and buyer use their blockchain addresses as unique IDs to join the data exchange and payment activities. The cryptocurrency based micro-payments protocol is designed to collect data in exchange for payment or other incentives. The key transaction records, such as invoices and receipts, are stored in a distributed ledger, so that all participants could verify the transactions.

    \noindent \emph{2. Identity Verification}: Since each blockchain account is uniquely indexed by its address that is derived from his/her own public key, the DDM administrator could deploy smart contract enabled microservices, which expose RESTful API to other microservices-enabled providers for referring identity verification results.

    \noindent \emph{3. Participants Rating}: To evaluate the quality of the data product provided by a given seller and reliability of payment by a buyer, rating smart contract based microservices are developed and deployed by trusted agents of the DDM system to allow buyers and sellers to rate each other. Owing to data transparency and immutability in the distributed ledger, the rating records stored in the blockchain network cannot be falsified or tampered with. The rating information would be used as recommended references for future potential buyers and sellers.

    \noindent \emph{4. Data Integrity}: Data integrity technologies are mainly to ensure data access at the same time avoid storing a huge amount of data in the blockchain. The data integrity microservices provides the dynamic data synchronization and efficient verification through a hashed index authentication process. The data publishers just simply interact with authorized ABI functions of the smart contract to save the hashed index of data to blockchain. In the hashed index verification process, data subscribers just query a hashed key-value index from the blockchain and compares it with calculated hash values of the received data.

    \noindent \emph{5. Access Control}: The data owners could transcode access control models and policies into a smart contract-based access control (AC) microservice. The AC microservices allows those service providers to control their data and resource without relying on a centralized third party authority, and it provides a decentralized access control solution to the DDM system.

    \noindent \emph{6. Privacy Policies}: The privacy microservices is mainly for privacy-sensitive data management, so that the privacy-sensitive data is not linked or accessed by unauthorized entities. Secure communications using encryption and saving sensitive data as a hash value mapping to actual data are efficient solutions to protect privacy. Furthermore, properly designing access control policies could effectively prevent unauthorized access to sensitive data.

\section{Conclusions}
\label{sec:conclusion}
In this paper, BlendSM-DDM, a blockchain-enabled decentralized security microservices framework is proposed to ensure data exchanges and payments in DDM systems like the Intelligent IoT Integrator (I3). The comprehensive overview of the system architecture is introduced and key elements are illustrated. The BlendSM-DDM could offer a scalable, flexible and lightweight security and privacy protection mechanism for business operations in DDM system. In ongoing work, we are designing the various elements presented in the architecture, and implementing a prototype running on current DDM system to verify the feasibility of the proposed solutions.

\bibliographystyle{IEEEtranS}