\section{Introduction}
As the development of modern deep learning techniques, more and more cutting-edge neural network structures have been widely used in various research and application fields. These machine learning frameworks in return inspire researchers from different background to design better neural network architectures to improve the model performance and solve more complex scenarios. One of the promising applications, quantitative finance, has recently been taking advantage of AI techniques to make lots of innovative algorithmic trading strategies for financial investment and portfolio optimization.

 Among all the machine learning frameworks, reinforcement learning (RL) has an unique advantage that the interactive training process is aligned with human decision making. During the past years, two basic optimization ways, Q-learning and Policy gradient have been developed to solve many RL problems, as well as the combination of the two. Take financial applications for example, Pastore et al.~\cite{2016Modelling} analyze data on 46 players from online games in financial markets and test whether the Q-Learning could capture these players'  behaviours based on a riskiness measure. Their results indicate that not all players are short-sighted, which contradicts with the naive-investor hypothesis. Besides, Li et al. (2019)~\cite{2019An} study the benefits of the three different classical deep RL models DQN~\cite{mnih2013playing}, Double DQN~\cite{van2016deep} and Dueling DQN~\cite{wang2016dueling} by predicting the stock price and conclude that DQN has the best performance. In addition, some researchers simulate and compare the improved deep RL method with the Adaboost algorithm and propose a hybrid solution. Interestingly, Lee et al. (2019)~\cite{2019Global} apply CNN to a deep Q-network which takes stock prices and volumes chart images as input to predict global stock markets. On the policy gradient side, Kang et al. (2018)~\cite{2018An} utilize the state-of-art Asynchronous Advantage Actor-Critic algorithm (A3C~\cite{mnih2016asynchronous}) to solve the portfolio management problem and design a standalone deep RL model. Moreover, Li et al. (2019)~\cite{2019Optimistic} propose an adaptive deep deterministic RL method (Adaptive DDPG) for some portfolio allocation tasks, which incorporates optimistic or pessimistic deep RL algorithm based on the influence of prediction error. Through analyzing data from Dow 30 component stocks, they show that the trading strategy outperforms traditional DDPG method~\cite{lillicrap2015continuous}. Sutta (2019)~\cite{2019Robust} compares the performance of an AI agent to the results of the buy-and-hold strategy and the expert trader by testing 15 years' forex data market with a paired t-Test. The findings show that AI can beat the buy \& hold strategy and commodity trading advisor in FOREX for both EURUSD and USDJPY currency pairs. However, these algorithms are quite expensive to train, considering high risk of investment loss at early stage in real financial market environment, even using simulated domains (Mnih et al. (2015)~\cite{mnih2015human}, Schulman et al. (2017)~\cite{schulman2017proximal}). A promising direction for improving sample efficiency is to explore model-based reinforcement learning (MBRL) methods~\cite{2019Model}. Chua, Calandra et al. (2019)~\cite{chua2018deep} study how to bridge this gap by employing uncertainty-aware dynamical models. They propose a novel algorithm called PETS, which approaches the asymptotic performance of-several benchmark model-free algorithms, while requiring significantly fewer samples. Janner et al. (2019)~\cite{janner2019trust} propose model-based Policy Optimization (MBPO) algorithm to boost PETS through monotonic improvement at each step and get state-of-the-art performance.

However, reinforcement learning is still plagued with intricate and difficult issues in the financial world. First, unforeseen variables such as international conditions, national policies, and business news events have an impact on financial markets. Second, hidden interactions among different agencies can be unpredictably difficult to understand.
 As it is shown from lots of stock marketing analysis, agencies with large capital investment can sometimes cause dramatic price fluctuation, leading to marketing panic and instability among the public. Additionally, complex nonlinearity, sometimes even complex chaotic behaviors, as well as unknowns like non-Gaussian noise tend to be present in financial marketing data, which results indistribution shift~\cite{cai2020distributed} of financial time series data over time. Complex causal relationships exist between different stocks. Stavroglou (2019)~\cite{nf_stavroglou2019hidden} characterizes the interaction of the local spatiotemporal dynamics between the attractors of asset pairs ($M_X$ and $M_Y$) as pattern causality (PC).  In this study, we empirically examine using the PC technique if there is a meaningful causal relationship between experimental stocks from various marketplaces and whether the magnitude of that relationship affects. There are many efforts have been made by academics to address the processes behind the aforementioned issues. On one hand, various data denoise methods are utilized to preprocess the financial time series data. For example, Bao W, et al. (2017)~\cite{2017A} propose a new model for stock prediction, wavelet transforms for decomposing and eliminating noise of the stock price time series, which helps stacked auto-encoders (SAEs) for extracted high-level deep features and long-short term memory (LSTM) to better forecast the stock price and improve predictive accuracy and profitability performance. On the other hand, researchers also create many kinds of indexes to optimize portfolio gain. In short, all these make it hard for us to build an reinforcement learning system which directly interacts with complex real marketing environment,  which motivates us to investigate model-based RL for algorithmic trading strategy optimization.

Despite the fact that stock prices are noisy and non-stationary, the difference between two consecutive days' stock prices follows a stationary distribution. In this article, we demonstrate statistically that the distribution of the alpha-stable levy will be followed by the difference. Classical MBRL algorithms, such as PETS and MBPO, often assume that each component of the state is independent of the others and follows a Gaussian distribution, which contradicts the complex hidden dependence and competition connections that occur in financial markets. Consequently, a model for the dynamic transfer process of the high-dimensional joint distribution of stocks must be developed using an algorithm.

By building a series of bijective transformations, normalizing flows (Dinh (2014)~\cite{dinh2014nice}; Dinh (2016)~\cite{dinh2016density}) proposed a mechanism for producing high-dimensional data distributions. Significant promise has already been shown in the creation and identification of images. Papamakarios et al.(2021)~\cite{papamakarios2021normalizing} describe flows using probabilistic modeling and inference. Kashif Rasul et al.(2020)~\cite{rasul2020multivariate} offer autoregressive deep learning architectures for simulating the temporal dynamics of multivariate time series using normalizing flows. Mevlana et al.(2016)~\cite{gemici2016normalizing } consider the problem of density estimation on Riemannian manifolds by employing normalizing flows.

All these excellent work inspires us to employ normalizing flows as the dynamic model of model-based RL algorithms for feature identification in complicated financial marketing environments, and to the best of our knowledge, rare literature has been found on utilizing MBRL into quantitative financial modeling.

In this paper,We simulate the alpha-stable levy distributions obeyed by the stock price difference series of several experimental stocks. Then, we employ normalizing flows for feature identification, and propose a novel strategy MBNF which is a novel model-based reinforcement learning algorithm by employing normalizing flows as dynamic model. Besides, we compare the performance of MBNF with MBPO method in three markets, and give the algorithm explanation by utilizing causal analysis and t-SNE method.

Overall, this paper makes the following major contributions:
\begin{itemize}
    \item\textbf{Marriage of NF and MBRL}.
    Normalizing flows is adopted to simulated high-dimensional joint probability of the complex trading environment, and we experimentally illustrate that normalizing flows will capture the relationships between stocks and make better decisions in the framework of model based reinforcement learning.

    \item\textbf{Business Performance}.
    We experiment various stocks from three different financial markets (Dow, NASDAQ and S\&P 500) with MBNF, and we show all these results of MBNF are comparatively better than modeling the state transition dynamics with independent Gaussian Processes.

    \item\textbf{Algorithm Explanation}.
    We utilize a causal analysis method to study the causal relationship among different stocks of the environment. Further, we print dynamic model loss of MBNF and MBPO  by visualizing high dimensional state transition data comparisons from real and virtual buffer with t-SNE, we uncover some effective patterns of better portfolio optimization strategies.
\end{itemize}

We structure our paper as the following sessions: First, in Section 2.1, we introduce some background knowledge of model-based RL, and then explain our financial environment setup under model-based RL framework in Section 2.2 and 2.3. We also interpret the transition dynamic from a dynamical system point of view and describe the Gaussian Process as some stochastic differential equation~(SDE), and correspondingly, the policy optimization of SAC~\cite{haarnoja2018soft} is connected with stochastic optimal control under the constraint of an SDE in Section 2.4. Next, we explain the normalizing flows method in detail in Section 3, and present the pseudo-code of our proposed model. Moreover, all the experiment results are shown in Section 4, as well as detailed explanations of how and why our proposed algorithm is better than classic MBPO. Finally, we summarize our conclusions and future research directions in Section 5.

\section{Framework}

Reinforcement learning aims to learn an optimal policy so as to maximize the expectation of cumulative rewards in the process of interacting with the environment.

\subsection{Space and Reward Definition}
By considering the investment trader as the intelligent agents and financial market as the corresponding environment, we model the stock trading problem as a Markov Decision Process (MDP). In this paper, we use bold to represent vectors

$\bullet \textbf{State space} \,\mathcal{S} = [B, \mathbf{P}, \mathbf{W}, \mathbf{I}]$:
State space $\mathcal{S}$ is the market information collected. $\forall s_t \in \mathcal{S}$ is the state of the agent at time $t$, which includes four parts information: account balance for agency $B_t$, current stock price $\mathbf{P_t}$, cumulative holding $\mathbf{W_t}$, technical indicators $\mathbf{I_t}$. Here the technical indicators $I_t$ is consisted of seven common technology factors: MACD~\cite{appel2003become}, SMA30, SMA60, BOLL~\cite{bollinger1992using}, RSI~\cite{ctuaran2011relative}, CCI~\cite{lai2020bidirectional}, ADX. Setting d stocks in each market, the table \ref{table_1} describes the notations of state $s_t$ in detail.

\begin{table}
\caption{Notations of State }
\label{table_1}
\centering
\begin{tabular}{l|lll}
\hline
Notation & Definition \\
\hline
$B_t$ & \begin{tabular}[c]{@{}l@{}} Account cash at time t;  $B_t \in \mathbb{R}_+$ \end{tabular} \\
\hline
$\mathbf{P_t}$ & \begin{tabular}[c]{@{}l@{}}  daily closing price at time t; $\mathbf{P_t} \in \mathbb{R}_+^{d}$  \end{tabular} \\
\hline
$\mathbf{W_t}$ & \begin{tabular}[c]{@{}l@{}} cumulative holding shares of each stock; $\mathbf{W_t} \in \mathbb{Z}_+^{d}$  \end{tabular} \\
\hline
MACD & \begin{tabular}[c]{@{}l@{}} Moving Average Convergence Divergence: \\a momentum indicator displays trend \end{tabular} \\
\hline
SMA30 & \begin{tabular}[c]{@{}l@{}} 30 day Simple Moving Average:\\ 30 day closing price equal weighted average   \end{tabular} \\
\hline
SMA60 & \begin{tabular}[c]{@{}l@{}} 60 day Simple Moving Average: \\60 day closing price equal weighted average   \end{tabular} \\
\hline
BOLL & \begin{tabular}[c]{@{}l@{}} Bollinger Bands:\\ judges the medium and long-term movement trend   \end{tabular} \\
\hline
RSI & \begin{tabular}[c]{@{}l@{}} Relative Strength Index: \\identifies inflection points of a trend  \end{tabular} \\
\hline
CCI & \begin{tabular}[c]{@{}l@{}} Commodity channel index: \\helps to find the degree of deviation of the price  \end{tabular} \\
\hline
ADX & \begin{tabular}[c]{@{}l@{}}  Average Directional Index:\\ determines the strength of a trend  \end{tabular} \\
\hline
\end{tabular}
\end{table}

$\bullet$ \textbf{Action Space} $\mathcal{A}$: Action space $\mathcal{A}$ is a set of available operations during the transaction on d stocks. The action taken by the agent at time $t$ denoted as $a_t$, is assumed to be finite dimensional and continuous. Here $a_t \in \mathcal{A}$ is a d-dimensional vector, where the $i_{th}$ dimension represents the action performed on the $i_{th}$ stock.

Denoting $\mathbf{W}^i_t$ as the cumulative shares of $i_{th}$ stock at time $t$. We assume the maximum trading volume for a single stock at one step is 100 shares. The available actions of the stock include buying, holding and selling:\\
• Buying: $a^i_t = + h$, $h$ shares can be bought and it leads to $\mathbf{W}^i_{t+1} = \mathbf{W}^i_t + h$.\\
• Selling: $a^i_t = - h$, $h$ shares can be sold from the current holdings. In this case, $\mathbf{W}^i_{t+1} = \max( \mathbf{W}^i_t - h, 0)$.\\
• Holding: $a_t = 0$, which means no change in $\mathbf{W}^i_t$.\\
where $h \in [0,100]$ is a positive integer.

$\bullet$ \textbf{Reward} $r_t$: The reward at time t is a mapping $ r_t: \mathcal{S} \times \mathcal{A} \rightarrow \mathbb{R}$. As is shown in equation (3), the direct reward $r_t$ of taking action $a_t$ in state $s_t$ is defined as the change amount of the net account value:

1)The total amount of the agent's net account value at time t, denoted by $Asset_t$, is the sum of current cash and the value of stock holdings. That is:
\begin{equation}
    Asset_t = B_t + \mathbf{P_t^T} \cdot \mathbf{W_t}
         = B_0 - \sum_{\tau = 0}^t \mathbf{P_{\tau}^T} \cdot a_{\tau} + \mathbf{P_{t+1}^T} \sum_{\tau = 0}^t a_{\tau}
\end{equation}

2) Considering the transaction cost at time t, we denote the cost as $C_t$:
\begin{equation}
    C_t = \mathbf{P_t^T} \cdot|\mathbf{a_t}|\cdot cost_{percentage}
\end{equation}
where $|\mathbf{a_t}|$ means the absolute value of each component of the vector $\mathbf{a_t}$.

3) We define the reward $r_t$ as a straightforward form:
\begin{equation}
    r_t = Asset_{t+1} - Asset_t - C_t
\end{equation}
where $B_t,\mathbf{P_t}$ is shown as Table 1; $\mathbf{a_t}$ is the action given by agent at time $t$; $\mathbf{P_t^T} \cdot \mathbf{W_t}$ is the inner product of vectors.

For RL problem, the agent’s objective is to maximize its expected reward. The reward can assess the quality of each action. What's more, the agent optimizes the strategy under the guidance of the reward. Noting the cumulative reward as $R_t$ in time t:
\begin{equation}
    R_t = \Sigma_{\tau=1}^t  r_{\tau}
\end{equation}

\subsection{Transition Dynamics}
We start with the definition of transition probability in the following.

$\bullet$ \textbf{Transition Probability} $\mathcal{P}$: The change of market conditions is considered as the state transition function. $\mathcal{P}: S \times A \times S \rightarrow [0,1]$ is a function of probabilities of state transitions. $$P_{s,s'}^a = \mathcal{P}(s_{t+1} = s' | s_t = s, a_t = a)$$.

In PETS, MBPO and M2AC, ensemble Gaussian Process is used to predict state transition dynamics. From stochastic analysis, we know that the state $s_t$ satisfies the following stochastic differential equation (SDE):
\begin{equation}
    dX_t = b(X_t,u_t)dt + \sigma(X_t,u_t)dW_t\label{SDE_control}
\end{equation}
where $b(X_t,u_t)$ and $ \sigma(X_t,u_t)$ are dirft term and difussion term. The initial condition $X_0 = x \in \mathbb{R}^d$, and $u_t \in \mathbb{R}^{d}$ is an $\mathbb{F}_t$-adapted control field and $\mathbf{W_t}$ is a $d$-dimensional $\mathbb{F}_t$-standard Brownian Motion.

$\bullet$ \textbf{Policy} $\pi$: Policy is a mapping characterized by the policy $\pi :\mathcal{S} \rightarrow \mathcal{A}$. $\forall t \in {1, . . . , T}$. The aim of policy $\pi$ is to maximize its total expected portfolio value.

The agent in state $s_t \in \mathcal{S}$ takes an action $a_t \in \mathcal{A}$ follows policy $\pi$, receives the reward $r_t = \mathcal{R}(s_t, a_t)$ and transits to the next state $s_{t+1}$ according to the transition probability $\mathcal{P}$. The state transition dynamics in PETS, MBPO, M2AC can be depicted as. The core optimization problem in RL is to find the optimal strategy $\pi^*$ that optimize total accumulated reward.

\subsection{Policy Optimization}
The actor-critic(AC) bridges the gap between value function approximation and policy gradient optimization for RL. It has been used for trading a large amount of stocks~\cite{xiong2018practical}. As a typical example of AC algorithms, the Deep deterministic policy gradient (DDPG) algorithm (Lillicrap et al. (2015)~\cite{lillicrap2015continuous}) shows great potential in stock trading problem and achieve some good results~\cite{emami2016deep,azhikodan2019stock}. However, it suffers from over sensitive to hyper-parameters. Tuomas Haarnoja proposed soft actor-critic (SAC) Algorithm~\cite{haarnoja2018soft} to solve the problem by employing entropy as the regularization term.

For random variable $x$ with probability density $p$, the entropy $H$ of $x$ can be defined as:
\begin{equation}
    H(p) = \underset{x\thicksim p}{\mathbb{E}}[-log p(x)]
\end{equation}
Considering discount factor $\gamma$, the entropy term is added to the expected reward, then the soft value function (soft Q function) becomes:
\begin{equation}
\begin{aligned}
    Q_{soft}^\pi(s, a)= & \underset{(s_t,a_t)\thicksim \rho_\pi}{\mathbb{E}} [\Sigma_{t = 0}^\infty \gamma^tr(s_t, a_t) + \\
    & \alpha \cdot \Sigma_{t=1}^\infty \gamma^t H(\pi(\cdot|s_t)) | s_0=s,a_0=a]
\end{aligned}
\end{equation}
where $\rho_\pi$ represents the distribution of state-action pair of the agent under policy $\pi$.

Taking equation (5) into account, we have the stochastic optimal control problem as follows:
\begin{equation}
    \mathcal{J}^u(x) = \mathbb{E}[\int_0^\infty f(X_s,u_s)e^{-\gamma s} ds |X_0^u = x]
\end{equation}
constrained on equation (\ref{SDE_control}).

Define
\begin{equation}
V(x) = \underset{u}{inf}\mathcal{J}^u(x) - \alpha log \pi(a|s)
\end{equation}

Then V satisfies the time-dependent HJB equation:
\begin{equation}
\underset{u}{inf}{\mathcal{L}^u V(x,u) + f(x,u) - \gamma V(x)} = 0
\end{equation}
where $\mathcal{L}^u V(x) = \frac{1}{2} Tr(\sigma \sigma^T Hess(V))(x,u) + b(x,u)^T \nabla V(x)$ is the generator of SDE (\ref{SDE_control}).

In Section 3.2, we will give the impact of normalizing flows on the HJB equations (10).

\section{Marriage of MBRL and Normalizing Flows}
\subsection{Curve Fitting}
The basic characteristics of our stock data are described in this subsection. Stock price differentials often follow a stationary distribution, despite the fact that stock prices are noisy and non-stationary. Based on the stationarity of price differences, the ARIMA model forecasts stock prices. Here we give the definition of smooth distribution.

An important property of the Gaussian distribution is that the sum of two independent and identically distributed random variables is still a Gaussian distribution.This property leads to the following definition.

\textbf{Definition 1.1} A random variable is called stable, if for any positive constant $a$ and $b$,  the following equation holds for some $c \in \mathbb{R}^{+}$ and $b \in \mathbb{R}$
$$
a X_1+b X_2 \stackrel{d}{=} cX+d
$$
Here $X_1$,$X_2$ are i.i.d with $X$,$\stackrel{d}{=}$ means equality in distribution.There is another equivalent definition of stable distribution.

\textbf{Definition 1.2} If for all $n \geq 2 $, there exist constants
$c_n>0$ and $d_n \in \mathbb{R}$ such that
$$
X_1+X_2+\cdots+X_n \stackrel{d}{=} c_n X+d_n
$$
where $X_i$ is i.i.d with $X$.For stable random variable, the characteristic function is in consistent format as following:
$$
\mathbb{E} \exp i \mu X= \begin{cases}\exp \left\{-\sigma^\alpha|\mu|^\alpha\left(1-i \beta(\operatorname{sign} \mu) \tan \frac{\pi \alpha}{2}\right)+i \gamma \mu\right\} & \text { if } \alpha \neq 1, \\ \exp \left\{-\sigma|\mu|\left(1+i \beta \frac{2}{\pi}(\operatorname{sign} \mu) \ln |\mu|\right)+i \gamma \mu\right\} & \text { if } \alpha=1 .\end{cases}
$$
where $\alpha \in(0,2], \beta \in[-1,1], \sigma \geq 0, \gamma \in \mathbb{R}$.
These four parameters determine a stable distribution function $S(\alpha, \sigma, \beta, \gamma)$ in which $\alpha$ is the tail index parameter indicating the fatness of the tail, $\beta$ the skewness parameter, $\gamma$ the scale parameter, and $\delta$ the location parameter.

$S(\alpha, \sigma, \beta, \gamma)$ differs in many respects from the common random variable.In most cases, $S(\alpha, \sigma, \beta, \gamma)$ can only be represented by a characteristic function except s.

We estimate the parameters of the difference between two consecutive days' stock prices on the training set (04/01/2011 to 04/07/2017). The data set presented here is a 1-day dataset; further details are provided in Section 4. For each stock, the difference between two consecutive days' stock prices $\delta \mathbf{P_t} = \mathbf{P_{t+\delta t}}-\mathbf{P_t}$, where $\delta t$ is the time interval.

The conventional MBPO model assumes that the state has a Gaussian distribution; however, if we consider fitting a stable distribution, we obtain different outcomes.
For consistency with our back-testing experiments data, without any processing of the differential data, Table \ref{DOW_fitting_appendix} \ref{SP_fitting_appendix} in Appendix A shows that the tail index parameter $\alpha$ are calculated broadly in between 1.5 and 1.75. The values are less than the $\alpha$=2 Gaussian distribution, indicating that difference prices are non-Gaussian with fat tails.

\begin{figure}[htbp]
    \begin{minipage}[b]{0.45\linewidth}
    \centering
    \includegraphics[width=1\textwidth]{NTES.pdf}
    \caption{NTES}
    \label{fit_NTES}
    \end{minipage}
    \begin{minipage}[b]{0.2\linewidth}
    \centering
    \begin{tabular}{ccccc}
    \hline
    \multicolumn{5}{c}{NAS}                           \\
    \hline
    ID & alpha  & beta    & mu    & sigma   \\
    \hline
    ASML       & 1.685  & -0.04   & 0.0627 & 0.7194  \\
    CMCSA      & 1.7411 & -0.0589 & 0.0206 & 0.1551  \\
    GOOGL      & 1.5915 & -0.0262 & 0.4345 & 3.8451  \\
    ISRG       & 1.6633 & 0.055   & 0.1117 & 1.462   \\
    MCHP       & 1.6925 & 0.0855  & 0.0233 & 0.3287  \\
    MSFT       & 1.5845 & 0.1665  & 0.0064 & 0.2581  \\
    NTES       & 1.2636 & 0.1     & 0.0108 & 0.1907  \\
    PEP        & 1.7476 & 0.045   & 0.0321 & 0.3594  \\
    TCOM       & 1.4954 & -0.0456 & 0.0123 & 0.3532  \\
    TMUS       & 1.6107 & -0.0379 & 0.0324 & 0.3287  \\
    \hline
    \end{tabular}
    \captionof{table}{NAS Market}
    \label{NAS_fitting}
    \end{minipage}%
\end{figure}

Figure \ref{fit_NTES} shows the fitted histogram using the stable distribution for NTES as example. The estimated stable distribution well characterizes the fat-tail behaviors and the bulk portion of difference price distributions.
It is worthy of mentioning that NTES appear to have $\alpha$ smaller than the other stocks, which is consistent with its fluctuation with prices skyrocketing and falling heavily at the beginning of 2017. The ITW have $\alpha$ close to 1.7, which is consistent with its fluctuation with prices skyrocketing and falling heavily at the beginning of 2017.

The four parameters $\alpha$, $\beta$, $\mu$ and $\sigma$, can offer clues to explain the main properties of difference price. However, $\beta$ and $\mu$ is not so robust to large price fluctuations and tends to have significant estimation errors. Still, the estimated $\beta$ and $\mu$ is close to 0, which means that difference price are not so skewed. The results provide additional views that price fluctuations for stock markets exhibit a symmetric behavior, which is also our finding. In this paper, we aim to improve the strategy returns by learning the combined levy distribution of high-dimensional differential pricing. We thus provide the NF model.

\subsection{Normalizing Flows}
Normalizing flows~\cite{dinh2016density} is a compound function which is consisted of a series of bijective transformations. It can generally be expressed as:
$$\mathbf{z} = F(\mathbf{x}) = f_n \circ f_{n-1} \circ f_{n-2} \circ \cdots \circ f_{2} \circ f_{1} (\mathbf{x})$$
where $\mathbf{z}$ follows a normal distribution whose density function is denoted as $p_{\mathbf{z}}(z)$ and $\mathbf{x}$ represents the distribution of the input data with density function $p_{\mathbf{x}}(x)$. $f_{i}$ is a simple invertible function.
then the sequence of random variable transformations can be written as a mapping flow as follows:
$$z_0=x\xrightarrow{f_1} z_1\xrightarrow{f_2}z_2\xrightarrow{f_3}\cdots\xrightarrow{f_n}z_n$$

Using the chain rule of derivation, we can express $p_{\mathbf{x}}(x)$ as
\begin{equation}
    p_{\mathbf{x}}(x)=p_{\mathbf{z}}(F(x))\left |  det \left(\frac{\partial F(x)}{\partial x}\right)\right |
\end{equation}
where $\frac{\partial F}{\partial x}$ is the Jacobi matrix.

An improved flow algorithm that combines addition and multiplication into the affining layer, sets the simple function $f_{i}$ as following:
$$
    \begin{cases}
        \mathbf{z}^{1:d}=\mathbf{x}^{1:d}\\
    \mathbf{z}^{d+1:D}=\mathbf{x}^{d+1:D}\odot exp(s(\mathbf{x}^{1:d}))+t(\mathbf{x^{1:d}})
    \end{cases}
$$

This simple function keeps the first d-dimensional variables unchanged, in which  $\odot$ is the element wise product, $s(\cdot)$ and $t(\cdot)$ are called scaling and translation function respectively which map $\mathbb{R}^d\rightarrow \mathbb{R}^{D-d}$.
In Real-NVP, the $s(\cdot)$ and $t(\cdot)$ of each affining layer are given by a neural network.

The Jacobi matrix of $f_{i}$ is
$$
\left(
    \begin{array}{cc}
        I_d & 0 \\
        \frac{\partial \mathbf{z}_2}{\partial \mathbf{x}_1} & diag(exp(s))
    \end{array}
\right).
$$

We can see that the determinant of this lower triangular matrix is $exp(\sum s)$, which means that the function is fully invertible. And the logarithm of the density function can be easily computed, as follows:

\begin{align}
   log\,p_{\mathbf{x}}(x)
   &= log\,p_{\mathbf{z}}(F(x))+\sum\limits_{i=1}^n log|det(\partial {f_i}/\partial x)|\nonumber \\
   &= log\,p_{\mathbf{z}} (F(x))+\sum\limits_{i=1}^n
   \sum\limits_{j=1}^d s_j(\mathbf{z}_{i-1}^{1:d})
\end{align}

The model is parameterized by scaling and translating the weights of the neural network $\theta$ and then trained on the training data points by stochastic gradient descent (SGD), where for each batch $D$, we maximize the average log likelihood (14) given by
\begin{align*}
    \mathcal{L}=\frac{1}{|\mathcal{D}|}\sum\limits_{\mathbf{x}\in\mathcal{D}}log\,p_x(\mathbf{x};\theta)
\end{align*}

\subsection{Proposed Model}
Based on normalizing flows as explained in the previous session, we know it is worthwhile to simulate high-dimensional joint probability of the complex trading environment. In \cite{lu2021learning}, the authors give a link between temporal normalizing folws of an SDE and and the solution of corresponding. For our case here, the density $\mathcal{P}_x(x)$ is the solution of Fokker-Plank equation corresponding to equation (5). This inspires us to simulate state transition dynamics through normalizing flows.

Our goal here is to predict next state $S_{t+1}$ from the current state $S_t$ and action $a_t$. In the prespective of optimal control, we discretize the transition dynamical system (5) in terms of Euler-Mariyam scheme:
\begin{equation}
X_{n+1} = X_n + b(X_n, u_n)\Delta t + \sigma (X_n, u_n)\sqrt{\Delta t}, \quad \quad n = 0,1,2,\cdots
\end{equation}
\begin{equation}
X_{n+1} - X_n = b(X_n, u_n)\Delta t + \sigma (X_n, u_n)\sqrt{\Delta t}, \quad \quad n = 0,1,2,\cdots
\end{equation}
If we sample from the probability density function learning from sequential data {$X_{n+1} - X_n$}, then by adding this sample with $X_n$, we could obtain a predicted state $X_{n+1}$.

From the last session, we know our main task is to simulate the probability density function $\mathcal{P}_X$ through normalizing flows by maximizing equation (14).
In fact, the $\mathcal{P}_X$ can help us to learn more information about equation (5). One way to study this information is through constructing an inverse problem and to learn unknown drift and diffusion terms in equation (5) with the method of maximum log-likelihood~\cite{fang2022end}. In this way, we could learn the intrinsic hidden dynamics from real data as a stochastic dynamical system, thus bridge a way to connect the relationship between optimal control and reinforcement learning. Once the unknown drift and diffusion coefficients are obtained, the generator $\mathcal{L}^u$ of equation (10) is also recovered so that the RL problem could be solved by solving the HJB equation.The pseudo-code is put in table \ref{table_19}.

\begin{table*}[!ht]
\renewcommand{\arraystretch}{1.3}
\caption{MBNF Framework}
\label{table_19}
\centering
\begin{tabular}{l}
\toprule
\textbf{Algorithm of MBNF : MBPO algorithm with Normalizing Flows(NF)} \\

\midrule
Initialize policy $\pi_\theta$, Normalizing Flows model $\Phi_{NF}$,  real environment buffer $\mathcal{B}_{env}$,\\
agent training buffer $\mathcal{B}_{agent}$,  real distribution $\mathcal{P}(x)$ \\

\bf{for 1:E steps do}\\
\quad  $\mathcal{B}_{env}$ goes one step further and add the data to itself\\
\quad \bf{for every M steps do} \\
\quad \quad Train Normalizing Flows model $\Phi_{NF}$ on $\mathcal{B}_{env}$ :\\
\quad \quad \bf{for 1:N steps do}\\
\quad \quad \quad $\Phi_{NF}$ computes approximate distribution $\widetilde{\mathcal{P}}$ through Maximum Likelihood Function: \\
\quad \quad \quad${\max}_{\widetilde{\mathcal{P}}(x)} \mathbb{E}_{x \sim \mathcal{P}(x)}   \lg \widetilde{\mathcal{P}}(x) $ \\
\quad \quad \quad where \\
\quad \quad \quad$\lg \widetilde{\mathcal{P}}(x) = -\frac{D}{2}\lg(2\pi) - \frac{1}{2}\| f(x) \|^{2} + \lg \vert \det[\frac{\partial f}{\partial x}] \vert$ with the bijective $f: \mathcal{P}(x) \to \Phi \sim \mathcal{N}(0,1)$ \\
\quad \quad Update agent training buffer $\mathcal{B}_{agent}$ :\\
\quad \quad Select a batch of data from $\mathcal{B}_{env}$ as current states $s_{t}$\\
\quad \quad Get a batch of action $a_t$ using policy $\pi_{\theta}$ \\
\quad \quad Get a batch of predicted states $\widetilde{s}_{t+1}$ and predicted $\widetilde{r}_{t+1}$ with trained $\Phi_{NF}$ using $a_t$ \\
\quad \quad Add $s_{t}$,$a_t$,$\widetilde{s}_{t+1}$,$\widetilde{r}_{t+1}$ to $\mathcal{B}_{agent}$ \\

\quad \bf{for 1:L steps do} \\
\quad \quad Get a batch of quadruplet $s_{t}$,$a_t$,$\widetilde{s}_{t+1}$,$\widetilde{r}_{t+1}$ from $\mathcal{B}_{agent}$ \\
\quad \quad Update policy parameters on $\mathcal{B}_{agent}$: $\theta\:\leftarrow\:\theta\:-\:\lambda_\pi\:\triangledown_\theta\:(\mathcal{J}_\pi(\theta, \mathcal{B}_{agent}))$\\
\bottomrule
\end{tabular}
\end{table*}

\section{Experiment}
In this session, we present some experimental results and corresponding explanations, from which we can see how and why normalizing flows helps improve the performance of classical model-based reinforcement learning model in a stock portfolio.

\subsection{Dataset}
All of our sample data are from Yahoo finance database.
To compare the effectiveness of our strategies across different markets, we conduct our experiments on three markets: the Dow Jones Industrial Average (Dow), National Association of Securities Dealers Automated Quotations (NASDAQ) and the Standard and Poor's 500 (S\&P 500). The Dow represents the trend of traditional U.S. industrial companies, while the NASDAQ favors technology companies and emerging sectors; The S\&P 500 is the most comprehensive of the three markets.

Considering all stocks available during the whole experiment period, the 10 stocks selected by their turnover rate within 180 days before Jan. 04, 2011 in the market respectively. The selected stocks have relatively low turnover rate so as to avoid big noise or fluctuations in unstable financial markets with small market capitalization. Here, as a liquidity indicator, the "turnover rate" refers to the frequency of stock changing hands in the market within a certain period of time. In this way, we can better explore the applicability and generalization capability of the proposed algorithm in this paper.

For each experiment, the dataset is historical daily price data from Jan. 4, 2011 to Mar. 31, 2021. Data from Jan. 4, 2011 to Jul. 3, 2017 (\textbf{1635 days}) is used as training set, and data from Jul. 4, 2017 to Jan. 30, 2019 (\textbf{396 days}) is used as validation set. The remaining data from Jan. 31, 2019  to Mar. 31, 2021 (\textbf{545 days}) is used as testing set. We train our agent on training data, then select model hyper-parameters by evaluating metrics such as annualized return, Sharpe ratio and maximum drawdown on validation data, and finally apply the selected model onto the testing data, which is shown in the following figures and tables. Notice that all the testing results are averaged from 10 random experiments.

Here we setup the following hyper-parameters for our experiments: initial balance $B_0 = 1e6$ dollars, and maximum shares per transaction $h_{max} = 100$.

\subsection{Portfolio Performance}
The performance of MBNF and MBPO strategies is compared in the three markets when the transition cost is ignored. All these strategies on testing data are compared with the Market Index baselines: DJIA, NDX and GSPC respectively.

The evaluation metrics, we use the following measurements to evaluate proposed model performance.
\begin{itemize}
    \item{Annualized Return: the geometric average amount of money earned by an investment strategy each year over a given time period.
    \item{Cumulative Return: reflects the overall effect of trading strategy in a certain time range.}
    \item{Annualized Volatility: annualized standard deviation of portfolio return shows the robustness of the agent.}
    \item{Sharpe Ratio~\cite{sharpe1998sharpe}: the return earned per unit volatility which is a widely-used measure of the performance of an investment.}
    \item{Calmar Ratio: the annualized return earned per unit maximum drawdown in the period.}
    \item{Stability: determines R-squared of a linear fit to the cumulative log returns.}
    \item{Maximum Drawdown: before a new peak is attained, the maximum loss from a peak-to-trough decline of an investment.}
\end{itemize}

\noindent\textbf{Overall performance with different measurements:}

In Table \ref{totle_return_table}, we use seven indicators to evaluate the performance of MBNF algorithms. It can be seen from the table that the proposed model MBNF exhibits better performance on the seven metrics in three markets, though the two strategies used for the experiment suffered huge disturbance from COVID-19. In the Dow market, MBNF and MBPO methods get the best performance. For example, MBNF improves baseline sharpe-ratio by twice as much as MBPO, and it outperforming baseline's annual return, sharpe-ratio, and stability metrics by 10\%, 0.34, and 0.64 respectively. The MBPO model outperforms the baseline's annual return, sharpe-ratio metrics are 4\% and 0.17. In S\&P 500 market, MBNF model has around 1.03 sharpe-ratio. In NASDAQ market, the performance of the two models is relatively poor, the baseline get the best score in the maximum drawdown metrics. The proposed MBNF model in annual return and sharpe ratio index is only slightly higher than the baseline. The MBPO model does not even meet the baseline benchmark. Besides, maximum drawdown metrics shows that the normalizing flows can help the agent to make a more stable investment decision with less riskiness.

And finally, considering the score of Dow and NASDAQ, we can see the S\&P 500 gets the middle situations in annual return, sharpe-ratio, stability, maximum-drawdown and other indicators, which also reflects that S\&P market stocks are more comprehensive.

Considering both risk management and profit gain, Sharpe ratio and Calmar ratio show that MBNF is significantly better than MBPO model and baseline, indicating combining normalizing flows with some classical reinforcement learning algorithms is a worthy direction to try.

\begin{table*}
\setlength{\tabcolsep}{3pt}
\centering
\caption{Performance comparison in three market}
\label{totle_return_table}
\scalebox{0.8}{
\begin{tabular}{c|ccc|ccc|ccc}
\hline
                      & \multicolumn{3}{c|}{Dow}                & \multicolumn{3}{c|}{NASDAQ}                   & \multicolumn{3}{c}{S\&P}                   \\
\hline
Backtest Indicators   & MBNF              & MBPO     & Baseline & MBNF              & MBPO     & Baseline          & MBNF              & MBPO     & Baseline  \\
\hline
Annualized Return     & \textbf{23.30\%}  & 17.58\%  & 14.00\%  & \textbf{35.35\%}  & 30.26\%  & 33.88\%           & \textbf{26.16\%}  & 23.93\%  & 19.48\%   \\
Cumulative Return     & \textbf{57.31\%}  & 41.95\%  & 32.69\%  & \textbf{92.43\%}  & 77.15\%  & 87.72\%           & \textbf{65.31\%}  & 59.04\%  & 46.85\%   \\
Annualized Volatility & \textbf{25.17\%}  & 24.02\%  & 26.59\%  & \textbf{28.66\%}  & 27.30\%  & 27.96\%           & \textbf{25.88\%}  & 25.22\%  & 25.25\%   \\
Sharpe Ratio          & \textbf{0.96}     & 0.80     & 0.63     & \textbf{1.20}     & 1.11     & 1.19              & \textbf{1.03}     & 0.98     & 0.83      \\
Calmar Ratio          & \textbf{78.32\%}  & 54.64\%  & 37.75\%  & \textbf{123.75\%} & 101.33\% & 120.84\%          & \textbf{91.52\%}  & 79.72\%  & 57.43\%   \\
Stability             & \textbf{90.33\%}  & 80.04\%  & 26.05\%  & \textbf{93.30\%}  & 85.59\%  & 89.53\%           & \textbf{83.22\%}  & 78.98\%  & 63.14\%   \\
Maximum Drawdown      & \textbf{-29.76\%} & -32.18\% & -37.09\% & -28.56\%          & -29.87\% & \textbf{-28.03\%} & \textbf{-28.59\%} & -30.01\% & -33.93\%  \\
\hline
\end{tabular}
}
\end{table*}

\noindent\textbf{Comparison of Daily Return:} To demonstrate the stability and effectiveness of our model, we show the mean, the maximum and the minimum daily return of ten experiments and the corresponding baseline in Figure~\ref{returns_dow}, Figure~\ref{returns_nas} and Figure~\ref{returns_sp}. We can see that the strategies showed a lot of volatility during COVID-19, the income of MBNF (purple) and its shadow shows robust advantages then MBPO (blue) and baseline (pink).

\begin{figure}[htbp]
	\centering
	\begin{minipage}{0.49\linewidth}
		\centering
		\includegraphics[width=0.9\linewidth]{returns_dow_true_0.pdf}
        \caption{Daily return curves for DOW.}
    \label{returns_dow}
	\end{minipage}
	\begin{minipage}{0.49\linewidth}
		\centering
		\includegraphics[width=0.9\linewidth]{returns_nas_true_0.pdf}
        \caption{Daily return curves for NAS.}
    \label{returns_nas}
	\end{minipage}
\end{figure}

\begin{figure}[htbp]
	\centering
	\begin{minipage}{0.49\linewidth}
		\centering
		\includegraphics[width=0.9\linewidth]{returns_sp_true_0.pdf}
        \caption{Daily return curves for SP.}
    \label{returns_sp}
	\end{minipage}
	\begin{minipage}{0.49\linewidth}
		\centering
		\includegraphics[width=0.9\linewidth]{Account_value_full_cost_true_false_new.pdf}
        \caption{Daily return curves for SP without technical indicators.}
    \label{cost_true_false}
	\end{minipage}
\end{figure}

\noindent\textbf{Technical Indicators for Normalizing flows:}
In the previous experiments, we add various technical factors into state space to better simulate the market environment, but normalizing flows can capture the relationships between stocks. Therefore, in order to further analyze the robustness of our method, we take the transition cost into account and test the performance of MBNF and MBPO with and without technical indicators in the state space of the S\&P market. As shown in Figure \ref{cost_true_false}, the performance of MBNF with and without the technical factor is not very different, but the MBPO model differs significantly better when the technical factor is used than when it is not. That means normalizing flows can discover the relationship between input data and make good decisions with it, but the gaussian model is more dependent on factors for trading due to the assumption that the input dimensions are independent of each other.

\subsection{Performance analysis / Environment Analysis}
In this part, we examine the causal linkages between stocks and the dynamic model's behavior throughout the backtesting process, which quantitatively explain the experimental performance and robutness of the proposed model.

\subsubsection{Causual Analysis:}
Financial markets have long been analyzed as complex systems with asset pricing being the foundation upon which structures such as financial networks are studied \cite{Normalizing_flows_2010contagion}
\cite{nf_boginski2005statistical}. Knowing the type of asset price interactions is vital to planning effective investment strategies or micro-trading tactics \cite{nf_fabozzi2008portfolio}.

This paper uses the pattern causality (PC) method of Stavroglou (2019)~\cite{nf_stavroglou2019hidden} to verify that correlation exists between different stocks within the same market and this holds for all three markets.  The positive and negative correlations between stocks are shown as Appendix A.

Figure \ref{Causual_Analysis} show the positive and negative correlations between stocks with PC method. Sub-figure \ref{Positive relationships in markets} illustrates the positive correlation between stocks. From top to bottom, the three columns in Sub-figure \ref{Positive relationships in markets} respectively represent the Dow, NASDAQ and S\&P market, and the three rows from left to right respectively represent the stock correlation scores greater than 0.3, 0.5 and 0.7. The Dow market has the strongest correlation, the NASDAQ market is slightly lower, and the S\&P market is significantly lower than the first two (first from the bottom), as can be seen from Sub-figure \ref{Positive relationships in markets}. This also explains why the MBNF model does not outperform the MBPO model by much in the S\&P market. In the NASDAQ market, although the experimental results show that the return of the MBNF model is only slightly higher than the baseline, it has definitely surpassed that of the MBPO model, which, after all, does not even reach the baseline.

Sub-figure \ref{Negative relationships in markets} illustrates the negative correlations between stocks, taking an alignment consistent with Sub-figure \ref{Positive relationships in markets}. Among the three markets, the Dow market has the highest negative correlation, followed by the S\&P market and the NASDAQ market has the lowest. The analysis of the positive and negative correlations across markets validates the introduction of the normalizing flows model because of the interdependence and competition between the experimental stocks.

\subsubsection{Dynamic Model Loss Analysis:}
In order to verify the reliability of the algorithm, we analyze the convergence of the dynamic model. Since MBNF uses log-likelihood as loss, while MBPO model uses mse method to calculate loss, we illustrate the loss changes of the two models separately in Figure \ref{Dynamic Model Loss} . As shown in Sub-figure \ref{NF_model_loss}, loss with NF as the dynamic model showed stable convergence when there were 40 episodes. The final convergence of model loss is around 0 in Dow market, and around -350 in NASDAQ and S\&P markets.

As shown in Sub-figure \ref{Gaussians_model_loss}, the Gaussian dynamic model converged within 20 episodes in the three markets. However, the model loss of convergence is stable between 500 and 1000.

\subsubsection{Buffer Analysis:}
The prediction of $s_{t+1}$ will affect the decisions from agent. Therefore, we visualize high dimensional state transition data comparisons from $\mathcal{B}_{model}$ of MBNF and MBPO algorithms with the t-distributed Stochastic Neighborhood Embedding (t-SNE) method to explore the better patterns of portfolio.

The t-SNE method is a popular nonlinear visualization algorithm for dimensionality reduction by employing a manifold learning method. It originated in image processing and has been widely used in various fields.

Figure \ref{true buffer and virtual buffer analysis} is the true buffer and virtual buffer analysis with t-SNE method. As shown in Sub-figure \ref{predict_$S_{t+1}$}, we can see that $s_{t+1}$ generated by normalizing flows has wilder exploration than the Gaussian dynamic. Considering the distribution of delta of two models in Sub-figure \ref{predict_delta}, we find that the Gaussian dynamic model has a strong Gaussian distribution for delta and NF has a more mixed Gaussian distribution.

The pattern may instruct that broader sampling and more variation in distribution allows normalizing flows to produce better trading strategies.

\section{Conclusion}
In this paper, we combine normalizing flows, which is good at generating high-dimensional joint distribution, as a dynamic model with model-based reinforcement learning algorithm to conduct a portfolio back-testing experiment in the Dow, NASDAQ and S\&P markets. The result illustrates normalizing flows will capture the correlation between high-dimensional stock data, facilitating the decision-making process of many automated trading systems.
Besides, we used the causual analysis method to verify the correlation between stocks in different markets, which suggests that we use normalizing flows as dynamic model is reasonable. This also explains why, in the S\&P market, the proposed model has not gotton a significant improvement over the MBNF model. For the NASDAQ market, although the model in this paper get a similar performance with baseline in metrics, MBNF robustly outperforms the MBPO model.

Furthermore, the model loss with normalizing flows and Gaussian are analyzed to demonstrate the convergence of MBNF, and the high-dimensional state transition data of buffers are visualized to find some effective patterns for better portfolio optimization strategies.

However, there are still some challenges that are open to us. For example,  the environment from real financial markets is complex and hard to fully simulate, considering the non-stationarity of the stock data, and unpredictable hidden causal factors. Besides, as we know that there is a close relationship between RL and stochastic optimal control. After setting up this question as a stochastic optimal control problem, how to prove the solution’s uniqueness and properties in an appropriate functional space, is still unknown to us. Moreover, the analysis of the generalization error bound of model-based RL is also important for us to design more adaptive and stable RL algorithms.