Continuous-Time Path-Dependent Exploratory Mean-Variance Portfolio Construction

Zhou Fang Department of Mathematics, The University of Texas at Austin

(Feb 2023)

Abstract

In this paper, we present an extended exploratory continuous-time mean-variance framework for portfolio management. Our strategy involves a new clustering method based on simulated annealing, which allows for more practical asset selection. Additionally, we consider past wealth evolution when constructing the mean-variance portfolio. We found that our strategy effectively learns from the past and performs well in practice.

1 Introduction

The revolutionary work Selection (\APACyear1952) is considered as the beginning of modern portfolio management, which proposes a framework for constructing portfolios in single periods. Much research in portfolio management came out after then. Li \BBA Ng (\APACyear2000) considers the portfolio selection problem in multiple periods setting and Zhou \BBA Li (\APACyear2000) further studies the portfolio selection problem in a continuous-time setting, which makes portfolio management in high frequency possible. In Wang \BOthers. (\APACyear2018), an innovative exploratory continuous-time mean-variance framework is first introduced. It replaces the deterministic policy with a stochastic policy to adjust assets’ holdings, achieving the ideal trade-off between exploration and exploitation, and is considered more robust than the deterministic policy. This exploratory framework is further studied in a series of subsequent papers. Wang (\APACyear2019) extends the framework into a multi-assets setting. Jia \BBA Zhou (\APACyear2022\APACexlab\BCnt2), Jia \BBA Zhou (\APACyear2022\APACexlab\BCnt1), Jia \BBA Zhou (\APACyear2022\APACexlab\BCnt3) develops a more general framework for policy evaluation and policy improvement.

In this paper, we further extend the exploratory continuous-time mean-variance framework by first clustering assets into several groups and then constructing a mean-variance portfolio by taking past wealth processes into consideration. In other words, there are phases in our framework. The first one is the clustering phase, and the second is the portfolio construction phase. The motivation for adding a clustering phase is to exclude similar assets and thus reduce the complexity. It is not possible and also not necessarily to hold thousands of stocks at one time for small retail and institutional investors. In the second portfolio construction phase, the motivation for considering the past wealth process or past performance is we think past experience should be important for making a current decision. When making a current decision, one should base it on past lessons together with the current wealth level instead of basing it only on the current wealth level. One intuitive example will be if a person came into the stock market in mid-2020, that person probably would think making money in the stock market is easy and maybe start to use high leverage in expect to make a higher gain, therefore, this person probably had a huge loss in 2022. However, if a person came into the stock market in 2008, that person probably would become more cautious, and therefore be aware of the market environment such as the Feds’ interest rate policy.

Many asset clustering methods have been developed in the past decades. As pointed out in Tang \BOthers. (\APACyear2022), one necessary criterion that wasn’t in the previous literature is that correlations for two assets from the same cluster with any other asset should be similar. Following that paper’s spirit, we use the similarity metric proposed in Dolphin \BOthers. (\APACyear2021) instead of using correlation. It is because correlation sometimes can’t capture similarities between two assets, as pointed out in Dolphin \BOthers. (\APACyear2021) one asset can have negative returns, and the other asset can have high returns, but these two assets still are highly correlated. The similarity metric we use here takes cumulative returns into consideration and therefore can capture the pattern similarities between assets. The clustering method we use in this paper is a simulated-annealing-based method, which is inspired by Ludkovski \BOthers. (\APACyear2022). We define an energy function that depends on asset similarities within a group. To our best knowledge, it is the first literature that uses this new similarity metric to cluster assets by a simulated annealing clustering method. From the empirical results, this method performs well, and can even possibly help investors find assets that guarantee statistical arbitrage opportunities, which requires further studies.

There is a huge literature on portfolio management or portfolio optimization, but very few consider portfolio management problems in a path-dependent setting. The path-dependent setting in portfolio management means that the positions of assets should depend on the performance of past investments. Learning from the past is an important philosophy that shouldn’t be ignored. The path-dependent case has technical difficulties, and there were no suitable tools for a long time. Thanks to the groundbreaking paper Dupire (\APACyear2019) that introduces functional Ito calculus, which defines a new type of calculus on path space, and thus proposes a path version of the Ito formula, and the Feynmann-Kac formula. However, the path-dependent Hamilton-Jacobi-Bellman equation is very hard to solve numerically. Thanks to the recent development of deep learning, makes solving path-dependent Hamilton-Jacobi-Bellman equations possible. There are literature that uses neural networks to solve path-dependent PDE, such as Saporito \BBA Zhang (\APACyear2020), and Sabate-Vidales \BOthers. (\APACyear2020). In this paper, we use PDGM to solve the PDE numerically, which is a neural network architecture whose main components are LSTM and feedforward networks. LSTM, as one of the most classical neural networks, has the ability to process sequential information and is naturally suitable to summarize past history. The feed-forward neural can then model certain functional by using the summaries given by LSTM. It is worth pointing out that LSTM is designed to forget some information in the past. In the deep learning community, the transformer basically replaced the LSTM to process sequential information. We believe replacing the LSTM component by a transformer component can have better results in solving path-dependent PDE, which requires further studies.

2 Simulated Annealing Clustering

Subsection 2.1 introduces a new similarity metric, and subsection 2.2 introduces the clustering method based on the proposed similarity metric via simulated annealing.

2.1 Similarity Metric of Financial Time Series

As pointed out in Dolphin \BOthers. (\APACyear2021), two assets can be perfectly correlated, but their performances will be very different. It is because conventional correlation is a similarity m, correlation cannot be the only criterion to cluster assets with similar patterns assets, correlation cannot be the only criterion. In this paper, we use the metric proposed in Dolphin \BOthers. (\APACyear2021) to measure the similarity of two assets.

Assume there are two time-series $X=[x_{1}.x_{2},...,x_{n}]$ , and $Y=[y_{1},y_{2},...,y_{n}]$ , which represent returns of two assets, and the similarity between these two time-series is defined to be

\mathbf{sim}(X,Y)=\frac{w}{1+e(X,Y)}+(1-w)\tau(X,Y),\hskip 5.69046pt

(1)

e(X,Y)=\sqrt{\big{(}\underset{1\leq i\leq n}{\Pi}(1+x_{i})-\underset{1\leq i\leq n}{\Pi}(1+y_{i})\big{)}^{2}}

(2)

\tau(X,Y)=\frac{\underset{1\leq i\leq n}{\sum}x_{i}y_{i}}{\sqrt{{\underset{1\leq i\leq n}{\sum}}x_{i}^{2}}{\sqrt{\underset{1\leq i\leq n}{\sum}y_{i}^{2}}}}

(3)

In the above definition of similarity, $e(X,Y)$ measures the distance between the cumulative returns of two assets and $\tau(X,Y)$ is a modified version of correlation that measures the similarity of absolute performance instead of relative performance to mean returns. $w$ is a hyper-parameter to be determined. Empirically, the metric performs best when $w$ is between 0.4 and 0.6.

2.2 Simulated Annealing Clustering

The clustering procedure is based on simulated annealing. The key to the simulated annealing clustering method is the energy functions. The ideal energy function should control the size of each cluster, and let similar assets be in the same group. Let $C=\{C_{1},C_{2},...,C_{k}\}$ be a partition of assets. Inspired by Ludkovski \BOthers. (\APACyear2022), we propose a similar energy function

\displaystyle E(C)=\underset{C_{i}\in C}{\sum}\Big{(}1-\frac{\kappa}{|C|-1}\underset{X,Y\in C_{i}}{\sum}\mathbf{sim}(X,Y)\Big{)}

(4)

In the above definition, $|C|$ is the number of clusters, which should never be 1, otherwise, it is not a clustering method. $\kappa$ is a hyper-parameter to be determined, which put a soft constraint on the number of clusters.

The simulated annealing clustering method is as follows – initially, all assets are in the same group, and for each iteration, a perturbation operation is applied to move one asset from one cluster to another with some criterion of accepting the perturbation. The criterion of accepting perturbation is whether or not the energy function is decreased or not. If the energy function is decreased, then one should accept the perturbation. If the energy function is increased, then one should accept the perturbation under a certain probability related to the energy function.

The exact simulated annealing clustering method is as follows
Initialization: Choose an initial temperature: $T_{0}$ , final temperature: $T_{f}$ , cooling rate: $\alpha<1$ , and let all assets in the same cluster.
Step 1: Assume current step is $l^{\text{th}}$ step. Denote the current partition as $C^{l}$ , and apply perturbation operation to the current partition, which resulting a new partition, denote as $C^{l+1}_{new}$ . Notice that the current temperature is $T_{l}$
Step 2: $\Delta E^{l}=E(C^{l+1}_{new})-E(C^{l})$ .
If: $\Delta E^{l}<0$ , accept partition $C^{l+1}_{new}$ as the partition $C^{l+1}$ for next iteration
Else:

$C^{l+1}_{new}$ is accepted as $C^{l+1}$ with probability $\text{exp}(-\frac{\Delta E^{l}}{T_{l}})$

$C^{l}$ is accepted as $C^{l+1}$ with probability $1-\text{exp}(-\frac{\Delta E^{l}}{T_{l}})$

Step 3: Lowering the temperature based on cooling rate: $T_{l+1}=\alpha T_{l}$ . Step 4: If $T_{l+1}<T_{f}$ , end the clustering procedure. Otherwise, repeat steps 1 - 3

3 Exploratory Mean-Variance Framework

The exploratory mean-variance framework is first proposed in Wang \BOthers. (\APACyear2018) for one-dimensional asset dynamics and extends to multi-dimensional asset dynamics in Wang (\APACyear2019). We follow Wang (\APACyear2019) to give an introduction to the exploratory mean-variance framework in this section.

The portfolio management problem is forming a portfolio that minimizes the variance under expected returns.

	$\displaystyle\text{min }\mathbf{Var}(x_{t})$		(5)
	$\displaystyle\text{s.t. }\mathbb{E}[x_{t}]=z$		(6)

Assume the clustering phase is already completed, and there are $n$ clusters. We randomly select $n$ assets, by choosing one asset from each cluster, whose dynamics are denoted as $\{S_{t}^{1},S_{t}^{2},...,S_{t}^{n}\}$ . For each asset $i$ , assume that the underlying dynamics to be

dS_{t}^{i}=S_{t}^{i}(\mu_{t}^{i}dt+\underset{1\leq j\leq n}{\sum}\sigma_{t}^{ij}dW_{t}^{j})

Let $\boldsymbol{\mu}_{t}=(\mu_{t}^{1},\mu_{t}^{2},...,\mu_{t}^{n})$ , $\boldsymbol{\sigma}_{t}=(\sigma_{t}^{ij})$ , and $d\boldsymbol{W}_{t}=(dW_{t}^{1},dW_{t}^{2},...,dW_{t}^{n})$ denote the vector of drift rates, the covariance matrix, and the vector of n-dim independent Brownian motions respectively. Let $\mathbf{a_{t}}=(a_{t}^{1},a_{t}^{2},...,a_{t}^{n})$ denote the holdings (in dollar) of those $n$ assets at time $t$ . Therefore, the wealth process will be

dx_{t}^{\mathbf{a_{t}}}=(\boldsymbol{\mu}_{t}-r\boldsymbol{e}_{d})^{T}\boldsymbol{a}_{t}dt+\boldsymbol{a}_{t}^{T}\boldsymbol{\sigma}_{t}d\boldsymbol{W}_{t}

(7)

In this paper, the holdings/control is path-dependent, which means holdings/control $\boldsymbol{a}_{t}$ depend on the past history of the wealth process. To differentiate the current wealth level and the path of the wealth process, we use lowercase letters for the current wealth level and use uppercase letters for the path of the wealth process. For example, $x_{t}^{\boldsymbol{a}_{t}}$ means the wealth level at time $t$ under holdings/control $\boldsymbol{a}_{t}$ , while $X_{t}$ means a path of wealth process from time $0$ to time $t$ . Now, the policy $\boldsymbol{\pi}(\boldsymbol{a}_{t}|X_{t})$ is a probability density, meaning that holdings/control at time $t$ has probability density $\boldsymbol{\pi}(\boldsymbol{a}_{t}|X_{t})$ to be at $\boldsymbol{a}_{t}$ , based on a realized wealth process is $X_{t}$ .

Under the exploratory framework proposed in Wang \BOthers. (\APACyear2018), we have a new stochastic process that denotes the average performance of the wealth process is

{d\tilde{x}_{t}}^{\boldsymbol{\pi}}=(\boldsymbol{\mu}_{t}-r\boldsymbol{e}_{d})^{T}\boldsymbol{m}_{t}dt+\boldsymbol{m}_{t}^{T}\boldsymbol{\sigma}_{t}d\boldsymbol{W}_{t}+\sqrt{\mathbf{Tr}(\boldsymbol{\Sigma}_{t}^{T}\boldsymbol{C}_{t})}d\widetilde{W_{t}}

(8)

Here, we assume that policy $\boldsymbol{\pi}$ is appied at time $t$ . $\widetilde{W_{t}}$ is an independent Brownian motion, $\boldsymbol{\Sigma}_{t}=\boldsymbol{\sigma}_{t}\boldsymbol{\sigma}_{t}^{T}$ is a positive definite matrix. $\boldsymbol{m}_{t}$ , and $\boldsymbol{C}_{t}$ are as follows,

\boldsymbol{m}_{t}=\int_{\mathbb{R}^{n}}\boldsymbol{a}_{t}\boldsymbol{\pi}(\boldsymbol{a}_{t}|\widetilde{X_{t}})d\boldsymbol{a}_{t}

(9)

\boldsymbol{C}_{t}=\int_{\mathbb{R}^{n}}[\boldsymbol{a}_{t}-\boldsymbol{m}_{t}][\boldsymbol{a}_{t}-\boldsymbol{m}_{t}]^{T}\boldsymbol{\pi}(\boldsymbol{a}_{t}|\widetilde{X_{t}})d\boldsymbol{a}_{t}

(10)

The goal is to identify a policy $\boldsymbol{\pi}$ that minimizes the following objective function, where $w$ is the Lagrangian multiplier

\displaystyle\mathbb{E}\Big{[}(\tilde{x}_{T}^{\boldsymbol{\pi}}-w)^{2}+\gamma\int_{0}^{T}\int_{\mathbb{R}^{n}}\boldsymbol{\pi}(\boldsymbol{a}_{t}|\widetilde{X_{t}}^{\boldsymbol{\pi}})\log\boldsymbol{\pi}(\boldsymbol{a}_{t}|\widetilde{X_{t}}^{\boldsymbol{\pi}})d\boldsymbol{a}_{t}dt\Big{]}-(w-z)^{2}

(11)

This cost function is very similar to the cost function in Wang (\APACyear2019), but instead, the cost function is path-dependent. The superscript $\boldsymbol{\pi}$ indicates wealth processes are generated under the policy $\boldsymbol{\pi}$ . Now that we have path-dependent control and path-dependent cost function, we can’t simply apply classical stochastic control results here, a functional version of stochastic control is needed.

4 Path Dependent Stochastic Control

4.1 Functional Ito Calculus

Functional Ito calculus is a crucial and powerful tool developed in Dupire (\APACyear2019) to study path-dependent stochastic calculus problems. The following is a quick introduction to functional Ito calculus and some important consequences in the Dupire (\APACyear2019).

Let $\Lambda_{t}$ denote the set of cadlag paths up to time $t$ , and specifically, $\Lambda_{t}^{n}$ is the space of cadlag paths that takes values in $\mathbb{R}^{n}$ . $\Lambda_{t}^{n\times k}$ denotes the space of two cadlag paths, one takes values in $\mathbb{R}^{n}$ and the other takes value in $\mathbb{R}^{k}$ . In general, $\Lambda=\underset{t\in[0,T]}{\bigcup}\Lambda_{t}$ , and $f:\Lambda\to\mathbb{R}^{d}$ is a functional.

If $X_{t}$ is a path that takes value in $\mathbb{R}$ , and $X_{t}(u)$ indicates the value at time $u$ . We define vertical, and flat extensions as follows, (For the path that takes value in $\mathbb{R}^{n}$ , the above vertical and flat extensions are done on each dimension individually)

X_{t}^{h}(u)=\begin{cases}X_{u}&\text{if $0\leq u<t$}\\ X_{t}+h&\text{if $t=u$}\end{cases}

X_{t,\delta t}(u)=\begin{cases}X_{u}&\text{if $0\leq u<t$}\\ X_{t}&\text{if $t\leq u\leq t+\delta t$}\end{cases}

Let $f:\Lambda\to\mathbb{R}$ be a functional, then the partial derivative as follows, (for $f:\Lambda\to\mathbb{R}^{d}$ , the following definition is on each dimension individually)

\Delta_{t}f(X_{t})=\underset{\delta t\to 0^{+}}{\text{lim}}\frac{f(X_{t,\delta t})-f(X_{t})}{\delta t}

\Delta_{x}f(X_{t})=\underset{h\to 0}{\text{inf}}\frac{f(X_{t}^{h})-f(X_{t})}{h}

To understand the functional Ito formula, we need to introduce the metric on the path space $\Lambda$ . Let $X_{t}$ , and $Y_{s}$ be two paths in $\Lambda$ , and without loss of generality, assume $t\leq s$ . Their distance is defined to be

d_{\Lambda}(X_{t},Y_{s})=||X_{t,s-t}-Y_{s}||_{\infty}+(s-t)

(12)

After defining the metric, we can define the continuity of a functional.

A functional $f:\Lambda\to\mathbb{R}$ is $\Lambda$ -continuous at $X_{t}\in\Lambda$ if $\forall\epsilon>0$ , $\exists\delta>0$ , and $\forall Y_{s}$ such that $d_{\Lambda}(X_{t},Y_{s})<\delta$ , there is $|f(X_{t})-f(Y_{s})|<\epsilon$ . A functional is $\Lambda$ -continuous if it is $\Lambda$ -continuous at all paths.

A functional $f:\Lambda\to\mathbb{R}$ is in $\mathbb{C}^{1,2}$ if it is $\Lambda$ -continuous, $C^{2}$ in x, and $C^{1}$ in t, with its partial derivatives also be $\Lambda$ -continuous.

Theorem 4.1 (Functional Ito Formula)

Let $x_{t}$ be a continuous semi-martingale process, $f\in\mathbb{C}^{1,2}$ , and $X_{t}$ is a path of process $x_{t}$ . Then, for any $t\in[0,T]$ ,

\displaystyle f(X_{t})=f(X_{0})+\int_{0}^{t}\Delta_{t}f(X_{s})ds+\int_{0}^{t}\Delta_{x}f(X_{s})dx_{s}+\frac{1}{2}\int_{0}^{t}\Delta_{xx}f(X_{s})d\langle x\rangle_{s}

(13)

4.2 Functional Feynman-Kac Formula and Functional HJB equation

Now, assume that $X_{t}$ is a path for the wealth process. If we use policy $\boldsymbol{\pi}$ for the left time, the cost functional are as follow,

J(X_{t},\boldsymbol{\pi})=\mathbb{E}\Big{[}(\tilde{x}_{T}^{\boldsymbol{\pi}}-w)^{2}+\gamma\int_{t}^{T}\int_{\mathbb{R}^{n}}\boldsymbol{\pi}(\boldsymbol{a}_{s}|\widetilde{X_{s}}^{\boldsymbol{\pi}})\log\boldsymbol{\pi}(\boldsymbol{a}_{s}|\widetilde{X_{s}}^{\boldsymbol{\pi}})d\boldsymbol{a}_{s}ds\hskip 2.84544pt|\hskip 2.84544pt\widetilde{X_{t}}^{\boldsymbol{\pi}}=X_{t}\Big{]}-(w-z)^{2}

(14)

To simplify our computation, denote the following functional as

\widetilde{\boldsymbol{b}}(s,\widetilde{X_{s}},\boldsymbol{\pi})=(\boldsymbol{\mu}_{s}-r\boldsymbol{e}_{d})^{T}\boldsymbol{m}_{s}

(15)

\widetilde{\boldsymbol{\sigma}}(s,\widetilde{X_{s}},\boldsymbol{\pi})=\big{(}\boldsymbol{m}_{s}^{T}\boldsymbol{\sigma}_{s}\hskip 2.84544pt|\hskip 2.84544pt\sqrt{\mathbf{Tr}(\boldsymbol{\Sigma}_{s}^{T}\boldsymbol{C}_{s})}\big{)}

(16)

f(s,\widetilde{X_{s}},\boldsymbol{\pi})=\int_{\mathbb{R}^{n}}\boldsymbol{\pi}(\boldsymbol{a}_{s}|\widetilde{X_{s}}^{\boldsymbol{\pi}})\log\boldsymbol{\pi}(\boldsymbol{a}_{s}|\widetilde{X_{s}}^{\boldsymbol{\pi}})d\boldsymbol{a}_{s}

(17)

Thus, the dynamics, and the cost functional are

d\widetilde{x_{s}}^{\boldsymbol{\pi}}=\widetilde{\boldsymbol{b}}(s,\widetilde{X_{s}},\boldsymbol{\pi})dt+\widetilde{\boldsymbol{\sigma}}(s,\widetilde{X_{s}},\boldsymbol{\pi})(d\boldsymbol{W}_{t}\hskip 2.84544pt|\hskip 2.84544ptd\widetilde{W_{t}})^{T}

(18)

J(X_{t},\boldsymbol{\pi})=\mathbb{E}\Big{[}(\tilde{x}_{T}^{\boldsymbol{\pi}}-w)^{2}+\gamma\int_{t}^{T}f(s,\widetilde{X_{s}},\boldsymbol{\pi})ds\hskip 2.84544pt|\hskip 2.84544pt\widetilde{X_{t}}^{\boldsymbol{\pi}}=X_{t}\Big{]}-(w-z)^{2}

(19)

Let $\tau>t$ , and consider the following,

Y_{\tau}=J(\widetilde{X_{\tau}},\boldsymbol{\pi})+\gamma\int_{t}^{\tau}f(s,\widetilde{X_{s}},\boldsymbol{\pi})ds

(20)

Consider the following,

	$\displaystyle\mathbb{E}[Y_{\tau}\|X_{t}]$	$\displaystyle=\mathbb{E}\Big{[}\mathbb{E}[(\tilde{x}_{T}^{\pi}-w)^{2}+\gamma\int_{\tau}^{T}f(s,\widetilde{X_{s}},\boldsymbol{\pi})ds\|\widetilde{X_{\tau}}]+\int_{t}^{\tau}f(s,\widetilde{X_{s}},\boldsymbol{\pi})ds\Big{\|}X_{t}\Big{]}-(w-z)^{2}$
		$\displaystyle=\mathbb{E}\Big{[}\mathbb{E}\big{[}(\tilde{x}_{T}^{\pi}-w)^{2}+\gamma\int_{\tau}^{T}f(s,\widetilde{X_{s}},\boldsymbol{\pi})ds\big{\|}\widetilde{X_{\tau}}\big{]}\Big{\|}X_{t}\Big{]}+\mathbb{E}\Big{[}\int_{t}^{\tau}f(s,\widetilde{X_{s}},\boldsymbol{\pi})ds\Big{\|}X_{t}\Big{]}-(w-z)^{2}$
		$\displaystyle=\mathbb{E}\Big{[}(\tilde{x}_{T}^{\pi}-w)^{2}+\gamma\int_{\tau}^{T}f(s,\widetilde{X_{s}},\boldsymbol{\pi})ds\Big{\|}X_{t}\Big{]}+\mathbb{E}\Big{[}\int_{t}^{\tau}f(s,\widetilde{X_{s}},\boldsymbol{\pi})ds\Big{\|}X_{t}\Big{]}-(w-z)^{2}$
		$\displaystyle=\mathbb{E}\Big{[}(\tilde{x}_{T}^{\pi}-w)^{2}+\gamma\int_{t}^{T}f(s,\widetilde{X_{s}},\boldsymbol{\pi})ds\Big{\|}X_{t}\Big{]}-(w-z)^{2}$
		$\displaystyle=J(X_{t},\pi)$
		$\displaystyle=Y_{t}$

Therefore, we can see the functional $Y_{\tau}$ is martingale, by the functional Ito formula, the following can be derived

\widetilde{\boldsymbol{b}}(\tau,\widetilde{X_{\tau}},\boldsymbol{\pi})\Delta_{x}J(\widetilde{X_{\tau}},\boldsymbol{\pi})+\Delta_{t}J(\widetilde{X_{\tau}},\boldsymbol{\pi})+\frac{1}{2}\widetilde{\boldsymbol{\sigma}}(\tau,\widetilde{X_{\tau}},\boldsymbol{\pi})\widetilde{\boldsymbol{\sigma}}(\tau,\widetilde{X_{\tau}},\boldsymbol{\pi})^{T}\Delta_{xx}J(\widetilde{X_{\tau}},\boldsymbol{\pi})+\gamma f(\tau,\widetilde{X_{\tau}},\boldsymbol{\pi})=0

(21)

To simplify the notation, assume the policy $\boldsymbol{\pi}$ be applied at time $0$ , then at time $t$ , and there is a path $X_{t}=\widetilde{X}^{\boldsymbol{\pi}}_{t}$ . , and we have functional Feynman-Kac formula,

\widetilde{\boldsymbol{b}}(t,X_{t},\boldsymbol{\pi})\Delta_{x}J(X_{t},\boldsymbol{\pi})+\Delta_{t}J(X_{t},\boldsymbol{\pi})+\frac{1}{2}\widetilde{\boldsymbol{\sigma}}(t,X_{t},\boldsymbol{\pi})\widetilde{\boldsymbol{\sigma}}(t,X_{t},\boldsymbol{\pi})^{T}\Delta_{xx}J(X_{t},\boldsymbol{\pi})+\gamma f(t,X_{t},\boldsymbol{\pi})=0

(22)

Now, let $V(X_{t})=\underset{\boldsymbol{\pi}}{\text{inf}}\hskip 2.84544ptJ(X_{t},\boldsymbol{\pi})$ , the following is an informal derivation of path-dependent HJB equation. The rigorous derivation of path-dependent HJB equation for portfolio management problems faces lots of technical challenges, which is beyond the scope of this paper.

$\displaystyle V(X_{t})$	$\displaystyle=\underset{\pi}{\text{inf}}\hskip 2.84544ptJ(X_{t},\boldsymbol{\pi})$	(23)
	$\displaystyle=\underset{\pi}{\text{inf}}\hskip 2.84544pt\mathbb{E}\Big{[}(\tilde{x}_{T}^{\pi}-w)^{2}+\gamma\int_{t}^{T}f(s,\widetilde{X_{s}},\boldsymbol{\pi})ds\Big{\|}X_{t}\Big{]}-(w-z)^{2}$	(24)
	$\displaystyle=\underset{\pi}{\text{inf}}\hskip 2.84544pt\mathbb{E}\Big{[}\mathbb{E}\big{[}(\tilde{x}_{T}^{\pi}-w)^{2}+\gamma\int_{t+dt}^{T}f(s,\widetilde{X_{s}},\boldsymbol{\pi})ds\big{\|}\widetilde{X}_{t+dt}\big{]}\Big{\|}X_{t}\Big{]}-(w-z)^{2}+\mathbb{E}\Big{[}\gamma\int_{t}^{t+dt}f(s,\widetilde{X_{s}},\boldsymbol{\pi})ds\Big{\|}X_{t}\Big{]}$	(25)
	$\displaystyle=\underset{\pi}{\text{inf}}\hskip 2.84544pt\mathbb{E}\Big{[}(\tilde{x}_{T}^{\pi}-w)^{2}+\gamma\int_{t}^{T}f(s,\widetilde{X_{s}},\boldsymbol{\pi})ds\Big{\|}X_{t}\Big{]}-(w-z)^{2}$	(26)
	$\displaystyle=\underset{\pi}{\text{inf}}\hskip 2.84544pt\mathbb{E}\Big{[}\mathbb{E}\big{[}(\tilde{x}_{T}^{\pi}-w)^{2}+\gamma\int_{t+dt}^{T}f(s,\widetilde{X_{s}},\boldsymbol{\pi})ds\big{\|}\widetilde{X}_{t+dt}\big{]}-(w-z)^{2}\Big{\|}X_{t}\Big{]}+\mathbb{E}\Big{[}\gamma\int_{t}^{t+dt}f(s,\widetilde{X_{s}},\boldsymbol{\pi})ds\Big{\|}X_{t}\Big{]}$	(27)

Since the first term in the above equation depends on the policy applied from time $t+dt$ , which is independent from the policy applied during the interval $[t,t+dt]$ . In other words, the policy $\pi$ appears in the above equation is a combination of policy applied in the interval $[t,t+dt]$ , and the policy applied after time $t+dt$ . Therefore, the above equation becomes

V(X_{t})=\underset{\pi}{\text{inf}}\hskip 2.84544pt\mathbb{E}\big{[}V(X_{t+dt})\big{|}X_{t}\big{]}+\mathbb{E}\big{[}\gamma\int_{t}^{t+dt}f(s,\widetilde{X_{s}},\boldsymbol{\pi})ds\big{|}X_{t}\big{]}

(28)

Moving the left-hand side to the right, and letting $dt\to 0$ , the above equation becomes,

$\displaystyle 0$	$\displaystyle=\underset{\pi}{\text{inf}}\hskip 2.84544pt\mathbb{E}\big{[}V(X_{t+dt})-V(X_{t})\big{\|}X_{t}\big{]}+\mathbb{E}\big{[}\gamma\int_{t}^{t+dt}f(s,\widetilde{X_{s}},\boldsymbol{\pi})ds\big{\|}X_{t}\big{]}$	(29)
	$\displaystyle=\underset{\pi}{\text{inf}}\hskip 2.84544pt\mathbb{E}\big{[}\Delta_{t}V(X_{t})+\Delta_{x}V(X_{t})d\tilde{x}_{t}^{\boldsymbol{\pi}}+\frac{1}{2}\Delta_{xx}V(X_{t})d\langle\tilde{x}^{\boldsymbol{\pi}}\rangle_{t}\big{\|}X_{t}\big{]}+\gamma f(t,X_{t},\boldsymbol{\pi})dt$	(30)
	$\displaystyle=\underset{\pi}{\text{inf}}\hskip 2.84544pt\big{[}\Delta_{t}V(X_{t})+\widetilde{\boldsymbol{b}}(t,X_{t},\boldsymbol{\pi})\Delta_{x}V(X_{t})+\frac{1}{2}\widetilde{\boldsymbol{\sigma}}(t,X_{t},\boldsymbol{\pi})\widetilde{\boldsymbol{\sigma}}(t,X_{t},\boldsymbol{\pi})^{T}\Delta_{xx}V(X_{t})+\gamma f(t,X_{t},\boldsymbol{\pi})\big{]}dt$	(31)

Divide both sides by $dt$ , we have the following path-dependent HJB equation

\Delta_{t}V(X_{t})+\underset{\boldsymbol{\pi}}{\text{inf}}\hskip 2.84544pt\big{\{}\widetilde{\boldsymbol{b}}(t,X_{t},\boldsymbol{\pi})\Delta_{x}V(X_{t})+\frac{1}{2}\widetilde{\boldsymbol{\sigma}}(t,X_{t},\boldsymbol{\pi})\widetilde{\boldsymbol{\sigma}}(t,X_{t},\boldsymbol{\pi})^{T}\Delta_{xx}V(X_{t})+\gamma f(t,X_{t},\boldsymbol{\pi})\big{\}}=0

(32)

Consider the formula inside the $inf$ , plug in $\widetilde{\boldsymbol{b}}(t,X_{t},\boldsymbol{\pi})$ , $\widetilde{\boldsymbol{\sigma}}(t,X_{t},\boldsymbol{\pi})$ , and $f(t,X_{t},\boldsymbol{\pi})$ , we have the following formula,

\int_{\mathbb{R}^{n}}\Big{[}(\boldsymbol{\mu}_{t}-r\boldsymbol{e}_{d})^{T}\boldsymbol{a}_{t}\Delta_{x}V(X_{t})+\frac{1}{2}\boldsymbol{a}_{t}^{T}\boldsymbol{\sigma}_{t}\boldsymbol{\sigma}_{t}^{T}\boldsymbol{a}_{t}\Delta_{xx}V(X_{t})+\gamma\log\boldsymbol{\pi}(\boldsymbol{a}_{t}|X_{t})\Big{]}\boldsymbol{\pi}(\boldsymbol{a}_{t}|X_{t})d\boldsymbol{a}_{t}

(33)

To find a candidate for the optimal policy, the following equation has to be independent of $\boldsymbol{a}_{t}$ and only depend on $X_{t}$

L(\boldsymbol{a}_{t},X_{t},\boldsymbol{\pi})=(\boldsymbol{\mu}_{t}-r\boldsymbol{e}_{d})^{T}\boldsymbol{a}_{t}\Delta_{x}V(X_{t})+\frac{1}{2}\boldsymbol{a}_{t}^{T}\boldsymbol{\sigma}_{t}\boldsymbol{\sigma}_{t}^{T}\boldsymbol{a}_{t}\Delta_{xx}V(X_{t})+\gamma\log\boldsymbol{\pi}(\boldsymbol{a}_{t}|X_{t})

(34)

U(\boldsymbol{a}_{t},X_{t})=(\boldsymbol{\mu}_{t}-r\boldsymbol{e}_{d})^{T}\boldsymbol{a}_{t}\Delta_{x}V(X_{t})+\frac{1}{2}\boldsymbol{a}_{t}^{T}\boldsymbol{\sigma}_{t}\boldsymbol{\sigma}_{t}^{T}\boldsymbol{a}_{t}\Delta_{xx}V(X_{t})

(35)

Assume $L(\boldsymbol{a}_{t},X_{t},\boldsymbol{\pi})$ depend on the $\boldsymbol{a}_{t}$ , and by the fact that $L(\boldsymbol{a}_{t},X_{t},\boldsymbol{\pi})$ is continuous at every point $\boldsymbol{a}_{t}$ , then without loss of generality, there exist two regions $V_{1}$ , and $V_{2}$ such that $L(\boldsymbol{a}_{t},X_{t},\boldsymbol{\pi})-L(\boldsymbol{a}^{\prime}_{t},X_{t},\boldsymbol{\pi})>h$ for points $\boldsymbol{a}_{t}\in V_{1}$ , and $\boldsymbol{a}^{\prime}_{t}\in V_{2}$ . Without loss of generality, assume the volume of these two regions is equal to $V$ . Now, let $\epsilon$ be sufficiently small, and consider the upgraded $\boldsymbol{\pi}^{*}$ that $\boldsymbol{\pi}^{*}(\boldsymbol{a}_{t}|X_{t})=\boldsymbol{\pi}(\boldsymbol{a}_{t}|X_{t})+\epsilon$ for $\boldsymbol{a}_{t}\in V_{1}$ , and $\boldsymbol{\pi}^{*}(\boldsymbol{a}_{t}|X_{t})=\boldsymbol{\pi}(\boldsymbol{a}_{t}|X_{t})-\epsilon$ for $\boldsymbol{a}_{t}\in V_{2}$ , and $\boldsymbol{\pi}^{*}$ is same for the rest space. Therefore, the change of $L(\boldsymbol{a}_{t},X_{t},\boldsymbol{\pi})\boldsymbol{\pi}(\boldsymbol{a}_{t}|X_{t})$ for regions $V_{1}$ is

	$\displaystyle L(\boldsymbol{a}_{t},X_{t},\boldsymbol{\pi}^{})\boldsymbol{\pi}^{}(\boldsymbol{a}_{t}\|X_{t})-L(\boldsymbol{a}_{t},X_{t},\boldsymbol{\pi})\boldsymbol{\pi}(\boldsymbol{a}_{t}\|X_{t})$	(36)
$\displaystyle=$	$\displaystyle\epsilon U(\boldsymbol{a}_{t},X_{t})+\gamma(\boldsymbol{\pi}(\boldsymbol{a}_{t}\|X_{t})+\epsilon)\log(\boldsymbol{\pi}(\boldsymbol{a}_{t}\|X_{t})+\epsilon)-\gamma\boldsymbol{\pi}(\boldsymbol{a}_{t}\|X_{t})\log\boldsymbol{\pi}(\boldsymbol{a}_{t}\|X_{t})$	(37)
$\displaystyle=$	$\displaystyle\epsilon U(\boldsymbol{a}_{t},X_{t})+\gamma(\boldsymbol{\pi}(\boldsymbol{a}_{t}\|X_{t})+\epsilon)\big{[}\log\boldsymbol{\pi}(\boldsymbol{a}_{t}\|X_{t})+\log(1+\frac{\epsilon}{\boldsymbol{\pi}(\boldsymbol{a}_{t}\|X_{t})})\big{]}-\gamma\boldsymbol{\pi}(\boldsymbol{a}_{t}\|X_{t})\log\boldsymbol{\pi}(\boldsymbol{a}_{t}\|X_{t})$	(38)
$\displaystyle=$	$\displaystyle\epsilon U(\boldsymbol{a}_{t},X_{t})+\epsilon\gamma\log\boldsymbol{\pi}(\boldsymbol{a}_{t}\|X_{t})+\gamma\epsilon+\mathcal{O}(\epsilon^{2})$	(39)

Here, since the $\epsilon$ is sufficiently small, we use Taylor expansion for the $\log(1+\frac{\epsilon}{\boldsymbol{\pi}(\boldsymbol{a}_{t}|X_{t})})$ to the (37). Similarly, for region $V_{2}$ , the change of the $L(\boldsymbol{a}_{t},X_{t},\boldsymbol{\pi})\boldsymbol{\pi}(\boldsymbol{a}_{t}|X_{t})$ is

	$\displaystyle L(\boldsymbol{a}_{t},X_{t},\boldsymbol{\pi}^{})\boldsymbol{\pi}^{}(\boldsymbol{a}_{t}\|X_{t})-L(\boldsymbol{a}_{t},X_{t},\boldsymbol{\pi})\boldsymbol{\pi}(\boldsymbol{a}_{t}\|X_{t})$	(40)
$\displaystyle=$	$\displaystyle-\epsilon U(\boldsymbol{a}_{t},X_{t})+\gamma(\boldsymbol{\pi}(\boldsymbol{a}_{t}\|X_{t})-\epsilon)\log(\boldsymbol{\pi}(\boldsymbol{a}_{t}\|X_{t})-\epsilon)-\gamma\boldsymbol{\pi}(\boldsymbol{a}_{t}\|X_{t})\log\boldsymbol{\pi}(\boldsymbol{a}_{t}\|X_{t})$	(41)
$\displaystyle=$	$\displaystyle-\epsilon U(\boldsymbol{a}_{t},X_{t})+\gamma(\boldsymbol{\pi}(\boldsymbol{a}_{t}\|X_{t})-\epsilon)\big{[}\log\boldsymbol{\pi}(\boldsymbol{a}_{t}\|X_{t})+\log(1-\frac{\epsilon}{\boldsymbol{\pi}(\boldsymbol{a}_{t}\|X_{t})})\big{]}-\gamma\boldsymbol{\pi}(\boldsymbol{a}_{t}\|X_{t})\log\boldsymbol{\pi}(\boldsymbol{a}_{t}\|X_{t})$	(42)
$\displaystyle=$	$\displaystyle-\epsilon U(\boldsymbol{a}_{t},X_{t})-\epsilon\gamma\log\boldsymbol{\pi}(\boldsymbol{a}_{t}\|X_{t})-\gamma\epsilon+\mathcal{O}(\epsilon^{2})$	(43)

Therefore, the change of (31) is

	$\displaystyle\int_{V_{1}}\big{[}\epsilon U(\boldsymbol{a}_{t},X_{t})+\epsilon\gamma\log\boldsymbol{\pi}(\boldsymbol{a}_{t}\|X_{t})+\gamma\epsilon\big{]}d\boldsymbol{a}_{t}-\int_{V_{2}}\big{[}\epsilon U(\boldsymbol{a}_{t},X_{t})+\epsilon\gamma\log\boldsymbol{\pi}(\boldsymbol{a}_{t}\|X_{t})+\gamma\epsilon\big{]}d\boldsymbol{a}_{t}$	(44)
$\displaystyle=$	$\displaystyle\epsilon\int_{V_{1}}U(\boldsymbol{a}_{t},X_{t})+\gamma\log\boldsymbol{\pi}(\boldsymbol{a}_{t}\|X_{t})d\boldsymbol{a}_{t}-\epsilon\int_{V_{2}}U(\boldsymbol{a}_{t},X_{t})+\gamma\log\boldsymbol{\pi}(\boldsymbol{a}_{t}\|X_{t})d\boldsymbol{a}_{t}+\epsilon\gamma\int_{V_{1}}d\boldsymbol{a}_{t}-\epsilon\gamma\int_{V_{2}}d\boldsymbol{a}_{t}$	(45)
$\displaystyle=$	$\displaystyle\epsilon\int_{V_{1}}U(\boldsymbol{a}_{t},X_{t})+\gamma\log\boldsymbol{\pi}(\boldsymbol{a}_{t}\|X_{t})d\boldsymbol{a}_{t}-\epsilon\int_{V_{2}}U(\boldsymbol{a}_{t},X_{t})+\gamma\log\boldsymbol{\pi}(\boldsymbol{a}_{t}\|X_{t})d\boldsymbol{a}_{t}+V-V$	(46)
$\displaystyle=$	$\displaystyle\epsilon\int_{V_{1}}L(\boldsymbol{a}_{t},X_{t},\boldsymbol{\pi})d\boldsymbol{a}_{t}-\epsilon\int_{V_{2}}L(\boldsymbol{a}_{t},X_{t},\boldsymbol{\pi})d\boldsymbol{a}_{t}>\epsilon hV$	(47)

From the above derivation, one can see that if the policy depends on $\boldsymbol{a}_{t}$ , there is always an improved policy by slightly adjusting the probability density.

(\boldsymbol{\mu}_{t}-r\boldsymbol{e}_{d})^{T}\boldsymbol{a}_{t}\Delta_{x}V(X_{t})+\frac{1}{2}\boldsymbol{a}_{t}^{T}\boldsymbol{\sigma}_{t}\boldsymbol{\sigma}_{t}^{T}\boldsymbol{a}_{t}\Delta_{xx}V(X_{t})+\gamma\log\boldsymbol{\pi}(\boldsymbol{a}_{t}|X_{t})=C(X_{t})

(48)

Thus, the optimal policy is of the following form, where $A$ is constant

\boldsymbol{\pi}(\boldsymbol{a}_{t}|X_{t})=A(X_{t})\hskip 2.84544pt\text{exp}\Big{(}-\frac{1}{\gamma}\big{(}(\boldsymbol{\mu}_{t}-r\boldsymbol{e}_{d})^{T}\boldsymbol{a}_{t}\Delta_{x}V(X_{t})+\frac{1}{2}\boldsymbol{a}_{t}^{T}\boldsymbol{\sigma}_{t}\boldsymbol{\sigma}_{t}^{T}\boldsymbol{a}_{t}\Delta_{xx}V(X_{t})\big{)}\Big{)}

(49)

By the fact that $\boldsymbol{\pi}(\boldsymbol{a}_{t}|X_{t})$ is a probability distribution, we have

\boldsymbol{\pi}(\boldsymbol{a}_{t}|X_{t})=\frac{\text{exp}\Big{(}-\frac{1}{\gamma}\big{(}(\boldsymbol{\mu}_{t}-r\boldsymbol{e}_{d})^{T}\boldsymbol{a}_{t}\Delta_{x}V(X_{t})+\frac{1}{2}\boldsymbol{a}_{t}^{T}\boldsymbol{\sigma}_{t}\boldsymbol{\sigma}_{t}^{T}\boldsymbol{a}_{t}\Delta_{xx}V(X_{t})\big{)}\Big{)}}{\int_{\mathbb{R}^{n}}\text{exp}\Big{(}-\frac{1}{\gamma}\big{(}(\boldsymbol{\mu}_{t}-r\boldsymbol{e}_{d})^{T}\boldsymbol{a}_{t}\Delta_{x}V(X_{t})+\frac{1}{2}\boldsymbol{a}_{t}^{T}\boldsymbol{\sigma}_{t}\boldsymbol{\sigma}_{t}^{T}\boldsymbol{a}_{t}\Delta_{xx}V(X_{t})\big{)}\Big{)}d\boldsymbol{a}_{t}}

(50)

Thus, $\boldsymbol{\pi}(\boldsymbol{a}_{t}|X_{t})$ is Gaussian distribution, more specifically,

\boldsymbol{\pi}(\hskip 2.84544pt\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}\hskip 2.84544pt|X_{t})\sim\mathcal{N}\Big{(}-(\boldsymbol{\sigma}_{t}\boldsymbol{\sigma}_{t}^{T})^{-1}(\boldsymbol{\mu}_{t}-r\boldsymbol{e}_{d})\frac{\Delta_{x}V(X_{t})}{\Delta_{xx}V(X_{t})}\hskip 2.84544pt,\hskip 2.84544pt(\boldsymbol{\sigma}_{t}\boldsymbol{\sigma}_{t}^{T})^{-1}\frac{\gamma}{\Delta_{xx}V(X_{t})}\Big{)}

(51)

Plug the policy back into HJB equation, as shown in Wang (\APACyear2019), and the HJB equation becomes,

\Delta_{t}V(X_{t})-\frac{||\boldsymbol{\sigma}_{t}^{-1}(\boldsymbol{\mu}_{t}-r\boldsymbol{e}_{d})||^{2}}{2}\frac{[\Delta_{x}V(X_{t})]^{2}}{\Delta_{xx}V(X_{t})}+\frac{\gamma}{2}\bigg{[}d-d\hskip 2.84544pt\text{ln }\frac{2\pi e\gamma}{\Delta_{xx}V(X_{t})}+\text{ln }\text{det}(\boldsymbol{\sigma}_{t}^{T}\boldsymbol{\sigma}_{t})\bigg{]}=0

(52)

In path-dependent cases, it is impossible to write an analytic solution for equation $(28)$ . Solving equation $(28)$ numerically is the only way and should be sufficient for practitioners. It is natural to solve path-dependent PDE numerically using neural networks because the deep learning community has developed ways to deal with sequential information.

5 Deep Learning PDE Solver

To solve $(28)$ numerically, it is important to approximate the value function, $V(X_{t})$ , where $X_{t}$ is a path of wealth process up to time $t$ . However, there is no magic in deep learning. None of the deep learning algorithms can directly deal with continuous sample paths and output values. Therefore, we have to discretize the time and sample path, plug the discretized version into a neural network, and get the output value.

Perhaps the most famous neural network to capture long-term dependence is Long Short-Term Memory (LSTM), which is a special type of recurrent neural network but doesn’t successfully mitigate the problem of gradient exploding and gradient vanishing (Those two are the same thing). Plugging our discretized version of the sample path into LSTM, we can get outputs that represent the path information in a neat way. After capturing the path information, there is one more thing to do – modeling the value function based on those outputs that contain path information. It is an easier step, because of the universal approximation theorem, we can just use a feed-forward neural network. In the rest of this section, we follow the framework proposed in Saporito \BBA Zhang (\APACyear2020) to solve our equation.

5.1 Long Short-Term Memory

First, LSTM is not a specific neural network structure, but more like a building block instead. The structure of an LSTM block is as follows,

In the above block diagram, there are three inputs, $a_{i-1}$ , $c_{i-1}$ , and $x_{i}$ , and there are two outputs $a_{i}$ , and $c_{i}$ . Intuitively, one can think in this way, $a_{i-1}$ , as a single number, is a neat way to represent all the information for the first $i-1$ steps, so people call it ”long-term memory”. While for $c_{i-1}$ , the short-term memory, as the name indicated, recent several steps’ information has more influence on it. At last, $x_{i}$ represents the information one gets from the current step. Those three memories interact with each other to generate new long-term memory and short-term memory, $a_{i}$ , and $c_{i}$ .

Back to our equation, assume now we have a sample path of the wealth process $X_{t}$ , where $t\in[0,T]$ . We discretize the time interval into $N$ equal periods, $0=t_{0}<t_{1}<...<t_{N}=T$ , and $\delta t=t_{i+1}-t_{i}$ , which will be useful later. Let $x_{i}$ , where $i=0,1,2,...,N$ be the wealth at time $t_{i}$ . So, we have the following network that consist of $N+1$ blocks of LSTM. Here, $a_{-1}=0$ , and $c_{-1}=0$ represent there is no memory initially.

In summary, we discretize a sample path $X_{t}$ into $N$ pieces, and plug it into a neural network that consists of $N+1$ LSTM blocks, we get an output vector $\boldsymbol{a}=(a_{0},a_{1},...,a_{N})$ . The next thing to do is use the output vector to model the value function $V(X_{t})$ .

5.2 Value Function Approximation

To approximate the value function $v(X_{t})$ , where $t_{i}\leq t<t_{i+1}$ , we use the most well-known neural network — feed-forward neural network, or fully connected neural network with three inputs, $t_{i}$ , $x_{i}$ , and $a_{i-1}$ , and one output, denote as $\phi(t_{i},x_{i},a_{i-1};\theta^{f})$ , where $\theta^{f}$ means the parameters in the feed-forward network. The exact structure of the feed-forward network (the number of hidden layers, the width of each layer, and the exact activation function) will be discussed in the empirical studies, and won’t be important for the rest of the section.

To summarize the above procedure, denote $u(X_{t};\theta^{f},\theta^{lstm})=\phi(t_{i},x_{i},a_{i-1};\theta^{f})$ , where $\theta^{lstm}$ denotes parameters of LSTM neural networks. The partial derivatives of the value function can thus be represented as (here, we adapt the notation in Saporito \BBA Zhang (\APACyear2020))

	$\displaystyle u(X_{t_{i}}^{h};\theta^{f},\theta^{lstm})=\phi(t_{i},x_{i}+h,a_{i-1};\theta^{f})$		(53)
	$\displaystyle u(X_{t_{i},\delta t};\theta^{f},\theta^{lstm})=\phi(t_{i},x_{i},a_{i-1};\theta^{f})$		(54)
	$\displaystyle\Delta^{[\delta t]}_{t}u(X_{t_{i}};\theta^{f},\theta^{lstm})=\frac{u(X_{t_{i},\delta t};\theta^{f},\theta^{lstm})-u(X_{t_{i}};\theta^{f},\theta^{lstm})}{\delta t}$		(55)
	$\displaystyle\Delta^{[h]}_{x}u(X_{t_{i}};\theta^{f},\theta^{lstm})=\frac{u(X_{t_{i}}^{h};\theta^{f},\theta^{lstm})-u(X_{t_{i}};\theta^{f},\theta^{lstm})}{h}$		(56)
	$\displaystyle\Delta^{[h]}_{xx}u(X_{t_{i}};\theta^{f},\theta^{lstm})=\frac{u(X_{t_{i}}^{h};\theta^{f},\theta^{lstm})-2u(X_{t_{i}};\theta^{f},\theta^{lstm})+u(X_{t_{i}}^{-h};\theta^{f},\theta^{lstm})}{h^{2}}$		(57)

In this way, with the ”right” parameters $\theta^{f}$ , and $\theta^{lstm}$ , we can approximate $V(X_{t})$ by $u(X_{t};\theta^{f},\theta^{lstm})$ . In order to achieve the ”right” coefficients, we need to train the above model with simulated sample paths.

Consider the equation $(28)$ , to simplify the notations, let

	$\displaystyle A=\frac{\|\|\boldsymbol{\sigma}_{t}^{-1}(\boldsymbol{\mu}_{t}-r\boldsymbol{e}_{d})\|\|^{2}}{2}$		(58)
	$\displaystyle B=\frac{\gamma}{2}\big{[}d-d\text{\hskip 1.9919ptln\hskip 1.9919pt}2\pi e\gamma+\text{ln\hskip 1.9919pt}\text{det}(\boldsymbol{\sigma}_{t}^{T}\boldsymbol{\sigma}_{t})\big{]}$		(59)
	$\displaystyle C=\frac{\gamma d}{2}$		(60)

Thus, equation $(28)$ becomes

\Delta_{t}V(X_{t})-A\frac{[\Delta_{x}V(X_{t})]^{2}}{\Delta_{xx}V(X_{t})}+C\hskip 1.9919pt\text{ln}\hskip 1.9919pt\Delta_{xx}V(X_{t})+B=0

(61)

If there are $M$ simulated paths, and we discretize each sample path into $N$ pieces. Denote $X_{t}^{(j)}$ be the $j^{th}$ simulated path. To save the space, let $\theta=[\theta^{f},\theta^{lstm}]$ . then the loss function to minimize is

\displaystyle J_{N,M}(\theta)

\displaystyle=\frac{1}{M}\frac{1}{N}\sum_{j=1}^{M}\sum_{i=0}^{M}\big{(}\Delta_{t}^{[\delta t]}u(X_{t}^{(j)};\theta)-A\frac{[\Delta_{x}^{[h]}u(X_{t}^{(j)};\theta)]^{2}}{\Delta_{xx}^{[h]}u(X_{t}^{(j)};\theta)}+C\hskip 1.9919pt\text{ln}\hskip 1.9919pt\Delta_{xx}^{[h]}u(X_{t}^{(j)};\theta)+B\big{)}^{2}

(62)

Here, $\alpha$ is a hyper-parameter to be determined depending on the situation, the last term is in order to let the average of final values of the simulated paths be as close to the expected return as possible.

6 Empirical Studies

There are two phases of empirical studies, the first is the clustering phase, and the second is the portfolio construction phase.

In the clustering phase, the simulated annealing clustering method with $\kappa=0.0001$ , $w=0.5$ , $\alpha=0.99$ , $T_{0}=100$ , and $T_{f}=0.1$ is used to cluster stocks in the S&P 500 index based on 2020-2022 historic data into 25 groups. The procedure is repeated 100 times to minimize the energy function. Figures 1-4 show the cumulative returns of stocks from four different groups. The results indicate that most stocks in a group perform similarly, while there may be some exceptional stocks that perform differently. This finding suggests that the clustering method could be used for pair trading to identify statistical arbitrage opportunities if the size of each cluster is limited to 2-5 stocks. However, this research direction is not explored in this paper and left for future studies.

In the second phase, 25 stocks are randomly selected from each group and a mean-variance portfolio is constructed based on their historic data. Before constructing the portfolio, the drift rate and covariance matrix for each day must be estimated. In this paper, the estimation is based on the past 75 days of stock returns, which may not be the most accurate method. More accurate estimation can be achieved with high-frequency data and sufficient computation resources.

The parameters used in the trading strategy were set as follows: $\gamma=0.01$ , $d=25$ , and $r=0$ since the interest rate has been zero for most of the past two years. The initial wealth path was $[1.0,1.01]$ and the trading strategy began with 1 dollar. At the start of each period, the current wealth path $X_{t}$ was used to generate hundreds of simulated paths with the same drift rate and volatility. These simulated paths were then used to train the neural network, which produced the policy parameters for generating holdings for the next period. A new wealth point was added to the current wealth path after each iteration, therefore, a new neural network needs to be trained. Due to limited computational resources, the trading strategy was only tested for 30 consecutive trading days (30 iterations) at three different periods covering bull, bear, and volatile markets. Figure 5-7 shows that this trading strategy outperformed the S&P 500 index.

7 Conclusion

This paper introduces a novel asset clustering method and extends the exploratory mean-variance framework to the path-dependent case. To further improve this framework, possible enhancements include replacing LSTM with transformer to enhance the neural network’s structure, using high-frequency data to estimate model parameters, and adding constraints on leverage or portfolio asset percentages to enhance the framework’s robustness.

References

Dolphin \BOthers. (\APACyear2021) \APACinsertmetastardolphin2021measuring{APACrefauthors}Dolphin, R., Smyth, B., Xu, Y.\BCBL \BBA Dong, R. \APACrefYearMonthDay2021. \BBOQ\APACrefatitleMeasuring Financial Time Series Similarity with a View to Identifying Profitable Stock Market Opportunities Measuring financial time series similarity with a view to identifying profitable stock market opportunities.\BBCQ \BIn \APACrefbtitleInternational Conference on Case-Based Reasoning International conference on case-based reasoning (\BPGS 64–78). \PrintBackRefs\CurrentBib
Dupire (\APACyear2019) \APACinsertmetastardupire2019functional{APACrefauthors}Dupire, B. \APACrefYearMonthDay2019. \BBOQ\APACrefatitleFunctional itô calculus Functional itô calculus.\BBCQ \APACjournalVolNumPagesQuantitative Finance195721–729. \PrintBackRefs\CurrentBib
Jia \BBA Zhou (\APACyear2022\APACexlab\BCnt1) \APACinsertmetastarjia2022policy_b{APACrefauthors}Jia, Y.\BCBT \BBA Zhou, X\BPBIY. \APACrefYearMonthDay2022\BCnt1. \BBOQ\APACrefatitlePolicy evaluation and temporal-difference learning in continuous time and space: A martingale approach Policy evaluation and temporal-difference learning in continuous time and space: A martingale approach.\BBCQ \APACjournalVolNumPagesJournal of Machine Learning Research231541–55. \PrintBackRefs\CurrentBib
Jia \BBA Zhou (\APACyear2022\APACexlab\BCnt2) \APACinsertmetastarjia2022policy_a{APACrefauthors}Jia, Y.\BCBT \BBA Zhou, X\BPBIY. \APACrefYearMonthDay2022\BCnt2. \BBOQ\APACrefatitlePolicy gradient and actor-critic learning in continuous time and space: Theory and algorithms Policy gradient and actor-critic learning in continuous time and space: Theory and algorithms.\BBCQ \APACjournalVolNumPagesJournal of Machine Learning Research231541–55. \PrintBackRefs\CurrentBib
Jia \BBA Zhou (\APACyear2022\APACexlab\BCnt3) \APACinsertmetastarjia2022q{APACrefauthors}Jia, Y.\BCBT \BBA Zhou, X\BPBIY. \APACrefYearMonthDay2022\BCnt3. \BBOQ\APACrefatitleq-Learning in continuous time q-learning in continuous time.\BBCQ \APACjournalVolNumPagesarXiv preprint arXiv:2207.00713. \PrintBackRefs\CurrentBib
Li \BBA Ng (\APACyear2000) \APACinsertmetastarli2000optimal{APACrefauthors}Li, D.\BCBT \BBA Ng, W\BHBIL. \APACrefYearMonthDay2000. \BBOQ\APACrefatitleOptimal dynamic portfolio selection: Multiperiod mean-variance formulation Optimal dynamic portfolio selection: Multiperiod mean-variance formulation.\BBCQ \APACjournalVolNumPagesMathematical finance103387–406. \PrintBackRefs\CurrentBib
Ludkovski \BOthers. (\APACyear2022) \APACinsertmetastarludkovski2022large{APACrefauthors}Ludkovski, M., Swindle, G.\BCBL \BBA Grannan, E. \APACrefYearMonthDay2022. \BBOQ\APACrefatitleLarge Scale Probabilistic Simulation of Renewables Production Large scale probabilistic simulation of renewables production.\BBCQ \APACjournalVolNumPagesarXiv preprint arXiv:2205.04736. \PrintBackRefs\CurrentBib
Sabate-Vidales \BOthers. (\APACyear2020) \APACinsertmetastarsabate2020solving{APACrefauthors}Sabate-Vidales, M., Šiška, D.\BCBL \BBA Szpruch, L. \APACrefYearMonthDay2020. \BBOQ\APACrefatitleSolving path dependent PDEs with LSTM networks and path signatures Solving path dependent pdes with lstm networks and path signatures.\BBCQ \APACjournalVolNumPagesarXiv preprint arXiv:2011.10630. \PrintBackRefs\CurrentBib
Saporito \BBA Zhang (\APACyear2020) \APACinsertmetastarsaporito2020pdgm{APACrefauthors}Saporito, Y\BPBIF.\BCBT \BBA Zhang, Z. \APACrefYearMonthDay2020. \BBOQ\APACrefatitlePDGM: a neural network approach to solve path-dependent partial differential equations Pdgm: a neural network approach to solve path-dependent partial differential equations.\BBCQ \APACjournalVolNumPagesarXiv preprint arXiv:2003.02035. \PrintBackRefs\CurrentBib
Selection (\APACyear1952) \APACinsertmetastarselection1952harry{APACrefauthors}Selection, P. \APACrefYearMonthDay1952. \BBOQ\APACrefatitleHarry Markowitz Harry markowitz.\BBCQ \APACjournalVolNumPagesThe Journal of Finance7177–91. \PrintBackRefs\CurrentBib
Tang \BOthers. (\APACyear2022) \APACinsertmetastartang2022asset{APACrefauthors}Tang, W., Xu, X.\BCBL \BBA Zhou, X\BPBIY. \APACrefYearMonthDay2022. \BBOQ\APACrefatitleAsset selection via correlation blockmodel clustering Asset selection via correlation blockmodel clustering.\BBCQ \APACjournalVolNumPagesExpert Systems with Applications195116558. \PrintBackRefs\CurrentBib
Wang (\APACyear2019) \APACinsertmetastarwang2019large{APACrefauthors}Wang, H. \APACrefYearMonthDay2019. \BBOQ\APACrefatitleLarge scale continuous-time mean-variance portfolio allocation via reinforcement learning Large scale continuous-time mean-variance portfolio allocation via reinforcement learning.\BBCQ \APACjournalVolNumPagesarXiv preprint arXiv:1907.11718. \PrintBackRefs\CurrentBib
Wang \BOthers. (\APACyear2018) \APACinsertmetastarwang2018exploration{APACrefauthors}Wang, H., Zariphopoulou, T.\BCBL \BBA Zhou, X. \APACrefYearMonthDay2018. \BBOQ\APACrefatitleExploration versus exploitation in reinforcement learning: a stochastic control approach Exploration versus exploitation in reinforcement learning: a stochastic control approach.\BBCQ \APACjournalVolNumPagesarXiv preprint arXiv:1812.01552. \PrintBackRefs\CurrentBib
Zhou \BBA Li (\APACyear2000) \APACinsertmetastarzhou2000continuous{APACrefauthors}Zhou, X\BPBIY.\BCBT \BBA Li, D. \APACrefYearMonthDay2000. \BBOQ\APACrefatitleContinuous-time mean-variance portfolio selection: A stochastic LQ framework Continuous-time mean-variance portfolio selection: A stochastic lq framework.\BBCQ \APACjournalVolNumPagesApplied Mathematics and Optimization4219–33. \PrintBackRefs\CurrentBib

Refer to caption — Figure 1: cumulative returns of group 1 stocks

	$\displaystyle\mathbb{E}[Y_{\tau}\|X_{t}]$	$\displaystyle=\mathbb{E}\Big{[}\mathbb{E}[(\tilde{x}_{T}^{\pi}-w)^{2}+\gamma\int_{\tau}^{T}f(s,\widetilde{X_{s}},\boldsymbol{\pi})ds\|\widetilde{X_{\tau}}]+\int_{t}^{\tau}f(s,\widetilde{X_{s}},\boldsymbol{\pi})ds\Big{\|}X_{t}\Big{]}-(w-z)^{2}$
		$\displaystyle=\mathbb{E}\Big{[}\mathbb{E}\big{[}(\tilde{x}_{T}^{\pi}-w)^{2}+\gamma\int_{\tau}^{T}f(s,\widetilde{X_{s}},\boldsymbol{\pi})ds\big{\|}\widetilde{X_{\tau}}\big{]}\Big{\|}X_{t}\Big{]}+\mathbb{E}\Big{[}\int_{t}^{\tau}f(s,\widetilde{X_{s}},\boldsymbol{\pi})ds\Big{\|}X_{t}\Big{]}-(w-z)^{2}$
		$\displaystyle=\mathbb{E}\Big{[}(\tilde{x}_{T}^{\pi}-w)^{2}+\gamma\int_{\tau}^{T}f(s,\widetilde{X_{s}},\boldsymbol{\pi})ds\Big{\|}X_{t}\Big{]}+\mathbb{E}\Big{[}\int_{t}^{\tau}f(s,\widetilde{X_{s}},\boldsymbol{\pi})ds\Big{\|}X_{t}\Big{]}-(w-z)^{2}$
		$\displaystyle=\mathbb{E}\Big{[}(\tilde{x}_{T}^{\pi}-w)^{2}+\gamma\int_{t}^{T}f(s,\widetilde{X_{s}},\boldsymbol{\pi})ds\Big{\|}X_{t}\Big{]}-(w-z)^{2}$
		$\displaystyle=J(X_{t},\pi)$
		$\displaystyle=Y_{t}$