A Deep Neural Network Algorithm for Linear-Quadratic Portfolio
  Optimization with MGARCH and Small Transaction Costs
We analyze a fixed-point algorithm for reinforcement learning (RL) of optimal
portfolio mean-variance preferences in the setting of multivariate generalized
autoregressive conditional-heteroskedasticity (MGARCH) with a small penalty on
trading. A numerical solution is obtained using a neural network (NN)
architecture within a recursive RL loop. A fixed-point theorem proves that NN
approximation error has a big-oh bound that we can reduce by increasing the
number of NN parameters. The functional form of the trading penalty has a
parameter $\epsilon>0$ that controls the magnitude of transaction costs. When
$\epsilon$ is small, we can implement an NN algorithm based on the expansion of
the solution in powers of $\epsilon$. This expansion has a base term equal to a
myopic solution with an explicit form, and a first-order correction term that
we compute in the RL loop. Our expansion-based algorithm is stable, allows for
fast computation, and outputs a solution that shows positive testing
performance.
Andrew Papanicolaou, Hao Fu, Prashanth Krishnamurthy, Farshad Khorrami
A Deep Neural Network Algorithm for Linear-Quadratic Portfolio Optimization with MGARCH and Small Transaction Costs

###### Abstract

We analyze a fixed-point algorithm for reinforcement learning (RL) of optimal portfolio mean-variance preferences in the setting of multivariate generalized autoregressive conditional-heteroskedasticity (MGARCH) with a small penalty on trading. A numerical solution is obtained using a neural network (NN) architecture within a recursive RL loop. A fixed-point theorem proves that NN approximation error has a big-oh bound that we can reduce by increasing the number of NN parameters. The functional form of the trading penalty has a parameter \(\epsilon>0\) that controls the magnitude of transaction costs. When \(\epsilon\) is small, we can implement an NN algorithm based on the expansion of the solution in powers of \(\epsilon\). This expansion has a base term equal to a myopic solution with an explicit form, and a first-order correction term that we compute in the RL loop. Our expansion-based algorithm is stable, allows for fast computation, and outputs a solution that shows positive testing performance.

_Keywords--_ Hetereoskedasticity, MGARCH, Fixed-point algorithms, Reinforcement learning, Deep neural networks

###### Contents

* 1 Introduction
	* 1.1 Background and Review
	* 1.2 Results in this Paper
* 2 Model and Optimization Problem
	* 2.1 Two-Step Iteration Scheme
	* 2.2 Convergence to a Fixed-Point
	* 2.3 Neural Network Policy Approximation
* 3 Small-\(\epsilon\) Asymptotic Analysis & Implementation
	* 3.1 Expansion of \(\lambda_{t}\)
	* 3.2 Neural Network Algorithm
###### Abstract

We consider a \(\epsilon\)-term model for the problem of a linear-quadratic problem with a quadratic penalty on the order-\(\epsilon\) term. The problem is solved by a quadratic penalty on the order-\(\epsilon\) term.
### Background and Review

Other machine learning applications in finance include [15] where stochastic gradient descent (SGD) with deep NN architecture is used for computing prices of American options on large baskets of stocks, and in [16] where an RL approach is used to numerically solve high-dimensional backward stochastic differential equations related to finance. In [17], the authors utilized an LSTM network for predicting the price movement with daily S&P500 data1. The performance of the LSTM is mixed during different periods. The short-term trend prediction in the price movement on NASDAQ by the deep network was studied by [18]. The authors utilized features from both fundamental and technical analysis as the network input. [19] used a graph network to predict the stock price movement on S&P500 data. In [20], they demonstrate how adversarial learning methods can be used to automate trading in stocks. General methods from control theory have been applied for optimal trading decisions in [21, 22]. The effects of transaction costs and liquidity are well-studied ([1, 23, 24]). In particular, the "aim portfolio" description given in [2] has been a key result for the management of large funds. The discussed works are based on supervised learning. Therefore, the non-supervised learning approaches should also be studied as they may address more complicated problems.

Footnote 1: [https://en.wikipedia.org/wiki/S%26P_500#cite_note-history-17](https://en.wikipedia.org/wiki/S%26P_500#cite_note-history-17)

It was originally thought that combining NN with RL recursions would lead to instability ([11]) due to the approximation errors caused by using the neural networks. However, hypothetically, the accuracy of NN approximations can be improved to within arbitrarily close bounds (see [25, 26]). Additionally, the actor-critic approach has been shown in settings to have a stabilizing effect ([27, 28]). The actor-critic approach utilizes two independent deep networks to represent a policy function (the actor) that provides a group of possible actions for a given state and an evaluation function (the critic) that evaluates the action taken by the actor based on the current policy function. By alternatively updating the actor and critic with the given objective function, the two networks converge. The DDPG algorithm ([9]) takes the actor-critic approach for solving RL problems with continuous control space and is based on the DPG algorithm in [29]. The DQN method in [8] is an RL algorithm that utilizes and trains a deep network to represent the Q-value function. Then at each instant, the DQN takes the action based on the Q-value network. The DQN does not use an actor-critic approach. and can only address discrete action problems, whereas the DDPG, which extends from DQN, can also address continuous action problems. But the DDPG is more difficult to train and sometimes unstable as it uses an actor-critic approach. Reinforcement learning was utilized in many topics, such as wireless ([30]) and mobile edge computing networks ([31, 32]).

In portfolio management, [33] utilized the DDPG to optimize the cryptocurrency portfolio. [20] utilized both DDPG and proximal policy optimization (PPO) for portfolio management. [34] considers portfolio optimization for a transaction cost. However, the authors only considered a linear cost. In this work, we considered a quadratic transaction cost that is more practical. Similar to the problem we consider in this paper is the linear-quadratic regulator with uncertainty in the (constant) coefficients. There are provable bounds for identifying the minimum run-time required from an on-policy approximation, after which the observed data is used to refine the policy to within a given tolerance of the optimal ([35, 36]). This approach to uncertainty in on-policy learning is analyzed as an actor-critic approach in [11].

### Results in this Paper

The analyses in this paper show the effectiveness of RL and NN-based policy approximations for solving linear-quadratic programs with non-constant coefficients and small penalization on control. We consider an iterative scheme for the optimal controls, for which we can prove convergence to a fixed point under some reasonable conditions. We extend this fixed point argument to show that NN approximation errors will compound over time, but that their total will be of order big-oh in the magnitude of an error that is reduced as the number of NN parameters is increased. The problem is of practical interest because the MGARCH covariance process does not allow for an explicit solution, such as those seen in other linear-quadratic optimizations ([37, 38]).

For a faster implementation, we propose an algorithm that exploits the smallness of the transaction cost parameter. In particular, we write the solution as a series expansion in powers of the transaction costs parameter, which is small. Using only a few terms from this expansion we form a sub-optimal solution that tends toward the optimum as the parameter decreases toward zero. This approach is similar to a standard NN-based policy gradient but is more stable (i.e., the difference between our RL output and the ground-truth optimal output is bounded) because the series terms are fast and easy to compute, and because the NN approximation error is present only in the higher-order terms and therefore is an
order of magnitude smaller than it would've been without the expansion. In our experiments, we take the first two terms in this series: the based term that is equal to an explicitly computable myopic solution and a first-order correction term that we compute using an RL loop and NN functional approximation. We implement this approximation using a single network that contains only fully-connected layers without requiring dueling networks or special machinery, and in our studies on market data, we can see that the correction term provides some improvement compared to using only the myopic control, which is defined in (22). The key difference between our RL strategy and the myopic strategy is that the myopic strategy makes decisions without forecasting the future, whereas the RL strategy forecasts the future to make decisions. Overall, the contribution of our paper includes the following:

* We show the effectiveness of RL and NN-based policy approximations for solving linear-quadratic programs with non-constant coefficients and small penalization on control.
* We prove the convergence of an iterative scheme to a fixed point for the optimal controls under some reasonable conditions.
* We show that NN approximation errors will be of order big-oh bound.
* We propose an algorithm based on NN approximation that exploits the smallness of the transaction cost parameter for faster implementation.
* We evaluate our algorithm on both synthetic and historical market data, which shows positive testing results.

The remaining part of this paper is organized as follows: First, the problem is mathematically formulated. A solution with a two-step iteration scheme is proposed, and its convergence is analyzed. A practical method for implementing the solution using neural networks is proposed. A small-\(\epsilon\) analysis regarding the neural network solution is shown. Lastly, the method is evaluated on synthetic market data and historical market data, and the results are discussed.

## 2 Model and Optimization Problem

Let \(R_{t}\in\mathbb{R}^{n}\) denote the vector of \(n\)-many assets' returns realized at time \(t\), and let \(\Sigma_{t-1}\) be the covariance of \(R_{t}\) given the information immediately prior at time \(t-1\). An MGARCH model ([39, 40]) is the following:

\[R_{t+1} =\mu+Z_{t+1} \tag{1}\] \[\Sigma_{t+1} =CC^{\top}+A\Sigma_{t}A^{\top}+BZ_{t+1}Z_{t+1}^{\top}B^{\top} \tag{2}\]

where \(\mu\) is the conditional mean vector, \(Z_{t+1}\sim\text{iid}(0,\Sigma_{t})\), and where \(A\), \(B\) and \(C\) are \(n\times n\) matrices. The following condition ensures that \(\Sigma_{t}\) is always invertible:

**Condition 1**.: _The matrix \(C\) given in (2) is full rank so that \(CC^{\top}\) is invertible, that is, there is a positive lower bound \(c=\inf_{\|v\|=1}v^{\top}CC^{\top}v>0\)._

Asset prices are calculated by compounding the returns (1). For \(1\leq i\leq n\), the \(i^{th}\) asset's price is

\[S_{t+1}^{i}=h(S_{t}^{i},R_{t+1}^{i}) \tag{3}\]

where \(h\) is some known function. A typical choice of \(h\) is \(h(s,r)=s(1+r)\) as used in [39, 40], but for technical reasons, we will need to impose the following condition:

**Condition 2**.: _The function \(h\) in (3) is finitely bounded away from zero. That is, \(\|h\|_{\infty}=\sup_{s,r}|h(s,r)|<\infty\) and there exists constant \(\underline{s}>0\) such that \(\inf_{s,r}h(s,r)\geq\underline{s}\)._

Denote the \(\mathbb{R}^{n}\) vector of these prices as

\[\vec{S}_{t}=\begin{pmatrix}S_{t}^{1}\\ S_{t}^{2}\\ \vdots\\ S_{t}^{n}\end{pmatrix}\,.\]
Next, define the covariance matrix of the dollar returns,

\[P_{t}=\Psi_{t}\Sigma_{t}\Psi_{t}\]

where

\[\Psi_{t}=\text{diag}(\vec{S}_{t})\.\]

Let \(X_{t}\in\mathbb{R}^{n}\) denote a manager's holdings in assets (in contract units). The returns (in dollar units) on this portfolio are \(\sum_{i}(S_{t+1}^{i}-S_{t}^{i})X_{t}^{i}\), the expected value of these returns is \(\mu^{\top}S_{t}X_{t}\), and their variance is \(X_{t}^{\top}P_{t}X_{t}\). The portfolio manager has a control \(\{a_{t},t=1,2,3,\dots\}\) that she selects at time \(t\) to change \(X_{t}\). The manager's control should be optimal with respect to her mean-variance preferences,

\[V(x,s,p)=\sup_{a}\sum_{t=1}^{\infty}\delta^{t}\mathbb{E}\left[f(a _{t},X_{t},\vec{S}_{t},P_{t})\Big{|}X_{0}=x,\vec{S}_{1}=s,P_{1}=p\right]\] (4) s.t. \[f(a_{t},X_{t},\vec{S}_{t},P_{t})=\underbrace{\mu^{\top}\Psi_{t}X _{t}-\frac{\epsilon q(\vec{S}_{t},P_{t})}{2}a_{t}^{\top}\Psi_{t}a_{t}}_{\text {expected return}}\underbrace{\frac{\gamma}{2}X_{t}^{\top}P_{t}X_{t}}_{\text{ risk penalty}}\] \[X_{t}=X_{t-1}+a_{t}\]

where \(0\leq\delta<1\) is a discount factor, \(\epsilon>0\) is a (small) parameter, and the function \(q:\mathbb{R}^{n}\times\mathbb{R}^{n\times n}\rightarrow\mathbb{R}^{+}\) is a transaction cost (or liquidity penalty) that is higher at times when it is harder to trade and lower at times when there is plenty of liquidity. This form of transaction cost penalty was introduced in [1] and the relationship with volatility was shown in [3].

**Condition 3**.: _We assume \(q\) is bounded away from zero,_

\[1\leq q(s,p)\leq\mathcal{X}<\infty\qquad\forall(s,p)\in\mathbb{R}^{n}\times \mathbb{R}^{n\times n} \tag{5}\]

_where \(\mathcal{X}\) is a known constant._

As mentioned earlier, the penalty on trading should depend on the instantaneous value of the covariance matrix. Many works use principal component analysis (PCA) to study the relationship between degrees of freedom in the stock returns' covariance matrix and overall market volatility. For example, [6] observes that relatively few eigenvectors are needed to capture the majority of market variance during times of high market stress, thus resulting in wider bid/ask spreads and higher transaction costs; the relationship is reversed during times of low market stress. Based on this dynamic, we use the condition number to excite the transaction cost function when the market has reduced degrees of freedom. In the examples, we take \(q(s,p)=\text{cond}(\text{diag}^{-1}(s)p\text{diag}^{-1}(s))\), which is based on the empirical observation that losses in the S&P500 market index occur when there is a large spike in the condition number of the covariance matrix, as shown in Fig. 1. Additionally, we take \(\gamma=\frac{1}{W_{0}}\sum_{i}(\overline{\Sigma}^{-1}\mu)^{i}\) where \(\overline{\Sigma}\) is the covariance matrix estimated from market data and \(W_{0}\) is the value of the new capital, i.e., we set risk aversion so that the goal is for \(W_{0}\) to be invested in the market.

The problem formulated in (4) is similar to the mean-variance preferences problem in [2] but with the added non-constant cost from \(q(\sigma)\), where \(\sigma^{2}=\Sigma\) and \(\Sigma\) is the estimated covariance matrix. A similar type of financial control problem was considered in [41]. An effective way to analyze this system is to take a Hamiltonian approach and write it using a vector of Lagrange multipliers,

\[\sum_{t=1}^{\infty}\delta^{t}\mathbb{E}\Big{[}f(a_{t},X_{t},\vec{ S}_{t},P_{t})+\lambda_{t}^{\top}(X_{t}-X_{t-1}-a_{t})\Big{]}\] \[=\sum_{t=1}^{\infty}\delta^{t}\mathbb{E}\Big{[}f(a_{t},X_{t}, \vec{S}_{t},P_{t})-(\delta\lambda_{t+1}-\lambda_{t})^{\top}X_{t}-\lambda_{t}^ {\top}a_{t}\Big{]}+\lim_{t\rightarrow\infty}\mathbb{E}[\delta^{t+1}\lambda_{t +1}^{\top}X_{t}]-\lambda_{1}^{\top}X_{0}\]
where we have used the transversality condition ([42])

\[\lim_{t\rightarrow\infty}\mathbb{E}(\delta^{t+1}\lambda_{t+1}^{\top}X_{t})=0.\]

First-order conditions in \(a_{t}\) and in \(X_{t}\) yield a forward-backward system

\[X_{t} =X_{t-1}-\frac{1}{\epsilon q(\vec{S_{t}},P_{t})}\Psi_{t}^{-1} \lambda_{t} \tag{6}\] \[\lambda_{t} =\delta\mathbb{E}_{t}\lambda_{t+1}-\Psi_{t}\mu+\gamma P_{t}X_{t}\ ; \tag{7}\]

where \(\mathbb{E}_{t}\) denotes expectation conditional on the information observed up to time \(t\). However, in the real world, (6) and

Figure 1: X-axis: the date (month - day). Left Y-axis: condition number. Right Y-axis: log(SPY) where SPY denotes the S&P500 index ETF. Blue plot: condition number of the covariance matrix with time. Red plot: log of SPY with time.
(7) cannot be directly solved. Therefore, in this paper, we propose an iteration scheme for (6) and (7):

\[X_{t}^{(k+1)} =X_{t-1}^{(k+1)}-\frac{1}{eq(\vec{S}_{t},P_{t})}\Psi_{t}^{-1} \lambda_{t}^{(k)}\] \[\lambda_{t}^{(k+1)} =\left(I+\widetilde{P}_{t}\Psi_{t}^{-1}\right)^{-1}\left(\delta \mathbb{E}_{t}\lambda_{t+1}^{(k+1)}-\Psi_{t}\mu+\gamma P_{t}X_{t-1}^{(k+1)}\right) \tag{8}\]

where \(\widetilde{P}_{t}=\frac{\gamma}{\epsilon q(\vec{S}_{t},P_{t})}P_{t}\) and \(X_{0}^{(k)}=x_{0}\) for all rounds of iteration \(k\). Then, we implement RL using neural networks (NNs) to estimate the limiting fixed point from (8). The convergence of iterations in (8) depends on if the following condition holds:

**Condition 4**.: _There is a constant \(\Delta(\epsilon)<1\) such that_

\[\delta\left\|\left(I+\widetilde{P}_{t}\Psi_{t}^{-1}\right)^{-1}\right\|\leq \Delta(\epsilon)\;,\]

_for all t._

If Condition 4 holds, then we can prove that (8) converges to a unique fixed point. Condition 4 would always hold if \(P_{t}\) and \(\Psi_{t}\) commute, but this is an unrealistic condition. For the data used in this paper, we check empirically that in fact Condition 4 holds for \(\epsilon<1\). For theoretical purposes, because we are considering the small-\(\epsilon\) parameterization, the following proposition is useful for confirming Condition 4:

**Proposition 1**.: _Assume Condition 1, Condition 2 and Condition 3. If \(\epsilon\chi\|h\|_{\infty}<\gamma\underline{s}^{2}c\), then denoting \(\kappa=X\|h\|_{\infty}/(\gamma\underline{s}^{2}c)\), we have_

\[\sup_{t}\left\|\left(I+\widetilde{P}_{t}\Psi_{t}^{-1}\right)^{-1}\right\|\leq \frac{\epsilon\kappa}{1-\epsilon^{2}\kappa^{2}}\;, \tag{9}\]

_from which it follows that Condition 4 will hold for \(\epsilon\) small enough._

Proof.: (see Appendix). 

### Two-Step Iteration Scheme

The scheme in (8) is the basis for an iterative algorithm with two steps. The policy is given by \(\lambda^{(k)}\) and is used to generate \(X^{(k+1)}\). Then, upon observation of \(X^{(k+1)}\), the updated Lagrange multiplier \(\lambda^{(k+1)}\) is obtained via a fixed-point iteration,

\[\lambda_{t}^{(k^{\prime}+1)}=\left(I+\widetilde{P}_{t}\Psi_{t}^{-1}\right)^{- 1}\left(\delta\mathbb{E}_{t}\lambda_{t+1}^{(k^{\prime})}-\Psi_{t}\mu+\gamma P _{t}X_{t-1}^{(k+1)}\right)\;, \tag{10}\]

for fixed \(k\) and for \(k^{\prime}\rightarrow\infty\). Algorithmically, this can be described as a 2-step iterative procedure for finding a fixed point:

1. \(X_{t}^{(k+1)}=X_{t-1}^{(k+1)}-\Psi_{t}^{-1}\lambda_{t}^{(k)}/(\epsilon q( \vec{S}_{t},P_{t}))\),
2. \(\lambda_{t}^{(k+1)}=\lim_{k^{\prime}\rightarrow\infty}\lambda_{t}^{(k^{\prime})}\) where \(\lambda_{t}^{(k^{\prime})}\) is given by (10).

If we consider a finite-time version of (6) and (7) with terminal condition \(\lambda_{T+1}^{(k)}\equiv 0\) for all \(k\) (i.e., on \(X_{t}^{(k)}=X_{T}^{(k)}\) for all \(t>T\)), and given Condition 4, then we have a contraction mapping with \(\lambda_{t}^{(k^{\prime})}\) converging to a unique fixed point as \(k^{\prime}\rightarrow\infty\), namely, \(\lambda_{t}^{(k+1)}\) for all \(t\leq T\).

### Convergence to a Fixed-Point

The two-step iteration described in Section 2.1 looks for a fixed point of \(\lambda^{(k^{\prime})}\) for a given \(X^{(k+1)}\). In this section we present a theorem stating that the pair \((X^{(k)},\lambda^{(k)})\) given by (8) converges to a fixed point, thus confirming that the forward-backward system of (6) and (7) has a unique solution. We prove these results for the finite-time version of (6) and (7) with terminal condition \(\lambda_{T+1}^{(k)}\equiv 0\) for all \(k\).
**Theorem 1**.: _Consider the finite-time problem with terminal condition \(\lambda_{T+1}^{(k)}\equiv 0\) for all \(k\). If we assume Condition 1, Condition 2, Condition 3 and Condition 4, then the iterations of (8) will converge to a unique fixed point._

Proof.: (see Appendix). 

**Remark 1**.: _Accuracy of the approximation of the infinite-time problem by a finite-time problem can be proven with optimality bounds and a squeeze lemma as \(T\to\infty\)._

The main idea in the proof of Theorem 1 is to show a contraction in \(\mathbb{E}\|\lambda_{t}^{(k+1)}-\lambda_{t}^{(k)}\|\), thus confirming the existence of a unique fixed point. The initial step for setting up the proof is to write the following forward equation for the iteration difference,

\[\lambda_{t}^{(k+1)}-\lambda_{t}^{(k)}=\delta(I+\widetilde{P}_{t} \Psi_{t}^{-1})^{-1}\mathbb{E}_{t}\left[\lambda_{t+1}^{(k+1)}-\lambda_{t+1}^{(k) }\right]-\gamma(I+\widetilde{P}_{t}\Psi_{t}^{-1})^{-1}P_{t}\left(\sum_{t^{ \prime}=1}^{t-1}\frac{1}{\epsilon q(\vec{S}_{t^{\prime}},P_{t})}\Psi_{t^{ \prime}}^{-1}(\lambda_{t^{\prime}}^{(k)}-\lambda_{t^{\prime}}^{(k-1)})\right) \tag{11}\]

which is derived from (8) by differencing between iteration \(k+1\) and \(k\). The following quantity is important for convergence,

\[\Omega=\frac{\chi\|h\|_{\infty}}{\bar{\mathfrak{s}}}. \tag{12}\]

The quantity in (12) is important because from (11) we obtain the following inequality,

\[\underset{t\leq T}{\sup}\mathbb{E}\left\|\lambda_{t}^{(k+1)}- \lambda_{t}^{(k)}\right\|\leq\frac{\Omega}{1-\Delta(\epsilon)}\sup_{t\leq T} \sum_{t^{\prime}=1}^{t-1}\mathbb{E}\left\|\lambda_{t^{\prime}}^{(k)}-\lambda_ {t^{\prime}}^{(k-1)}\right\|\,\]

An equivalent form of this inequality appears in Theorem 1 and is used to show that \(\mathbb{E}\left\|\lambda_{t}^{(k+1)}-\lambda_{t}^{(k)}\right\|\to 0\) as \(k\to\infty\) for \(t\leq T\), thus proving the existence of a unique fixed point of (8) for finite \(T\).

### Neural Network Policy Approximation

The machine learning approach to solve (8) is to implement RL using neural networks (NNs). It amounts to estimating the limiting fixed point \(\lambda^{*}\) from (8) as

\[\lambda_{t}^{*}\approx\lambda(t,X_{t-1},\vec{S}_{t},P_{t};\theta)\]

where \(\lambda(\cdot,\cdot,\cdot,\cdot;\theta)\) is a policy approximation function with a feed-forward NN, and \(\theta\) denotes the NN parameters to be estimated. The NN takes into account time \(t\) so that the solution can be adapted to the time remaining until terminal time \(T\). An example of an NN that can accommodate time dependence is the Deep BSDE architecture in [16].

Given an initial estimate \(\theta^{(0)}\), we proceed to iteratively look for \(\theta\) that is close to a fixed point. The following iterative estimation scheme is the basis for the algorithm we'll implement:

\[\theta^{(k+1)} =\underset{\theta}{\arg\min}\ \ \sum_{t=1}^{T}\mathbb{E}\Big{\|} \lambda(t,X_{t-1}^{(k+1)},\vec{S}_{t},P_{t};\theta)-Y_{t}^{(k)}(\theta)\Big{\|} ^{2}\ \text{for}\] (13) s.t. \[X_{t}^{(k+1)} =X_{t-1}^{(k+1)}-\frac{1}{\epsilon q(\vec{S}_{t},P_{t})}\Psi_{t} ^{-1}\lambda(t,X_{t-1}^{(k+1)},\vec{S}_{t},P_{t};\theta^{(k)})\] \[Y_{t}^{(k)}(\theta) =\delta\left(I+\widetilde{P}_{t}\Psi_{t}^{-1}\right)^{-1}\mathbb{ E}_{t}\Big{[}\lambda(t+1,X_{t}^{(k+1)},\vec{S}_{t+1},P_{t+1};\theta)\Big{]} \mathbb{1}_{t<T}-\left(I+\widetilde{P}_{t}\Psi_{t}^{-1}\right)^{-1}\left(\Psi_ {t}\mu-\gamma P_{t}X_{t-1}^{(k+1)}\right)\]

where the indicator \(\mathbb{1}_{t<T}\) is used to enforce the finite-time problem's terminal condition of \(\lambda_{T+1}=0\). The following theorem proves the NN scheme in (13) will in fact be a decent approximation for the fixed point \(\lambda_{t}\).
**Theorem 2**.: _Assume Condition 1, Condition 2 and Condition 3. Furthermore, assume the family_

\[(\lambda(t,X_{t-1}^{(k)},\vec{S}_{t},P_{t};\theta^{(k)}))_{k=1,2,\ldots}\]

_are continuous and uniformly integrable, i.e., for any \(\eta>0\) there is compact set \(\mathbf{K}\subset\mathbb{R}^{d}\times\mathbb{R}^{d}\times\mathbb{R}^{d\times d}\) on which we have_

\[\sup_{k}\mathbb{E}\Big{[}\|\lambda(t,X_{t-1}^{(k)},\vec{S}_{t},P_{t};\theta^{( k)})\|\,\mathbb{1}_{\{(X_{t-1}^{(k)},\vec{S}_{t},P_{t})\notin\mathbf{K}\}}\Big{]} \leq\eta\;. \tag{14}\]

_Assume also that_

\[\mathbb{E}\Big{[}\|\lambda^{*}(t,X_{t-1}^{*},\vec{S}_{t},P_{t})\|\,\mathbb{1} _{\{(X_{t-1}^{*},\vec{S}_{t},P_{t})\notin\mathbf{K}\}}\Big{]}\leq\eta.\]

_For each \(k\) let \(\varepsilon^{(k)}\) denote the error from the NN approximation,_

\[\sup_{t\leq T}\!\!\mathbb{E}\Big{[}\,\Big{\|}\lambda(t,X_{t-1}^{(k)},\vec{S}_ {t},P_{t};\theta^{(k)})-Y_{t}^{(k-1)}(\theta^{(k)})\Big{\|}\,\mathbb{1}_{\{(X _{t-1}^{(k)},\vec{S}_{t},P_{t})\in\mathbf{K}\}}\Big{]}\leq\varepsilon^{(k)}\;.\]

_Then, the error of the iteration scheme in (13) is_

\[\sup_{t\leq T}\!\!\mathbb{E}\,\Big{\|}\lambda(t,X_{t-1}^{(k+1)},\vec{S}_{t},P_ {t};\theta^{(k+1)})-\lambda_{t}^{*}\Big{\|}=\mathcal{O}\left(\left(\sup_{\ell }\varepsilon^{(\ell)}+2\eta\right)\exp\Big{(}\frac{\Omega T}{1-\Delta(\epsilon )}\Big{)}\right)\;, \tag{15}\]

_for \(k\) large, where \(\lambda^{*}\) is the fixed point of (8)._

Proof.: (see Appendix). 

The bound in (15) is similar to the bounds derived in [43], wherein the approximation error was computed to be the sum of three terms: a sampling error term, an NN parameter error, and a value function estimation error. An additional similarity is an exponential growth in the bounding constants as \(T\) increases. In theory, the exponential growth in (15) is contained by continually increasing the hyperparameters so that \(\sup_{\ell}\varepsilon^{(\ell)}+2\eta\) tends to zero, but in practice increasing hyperparameters requires a growing number of training samples, leading to an infeasibly long computation time. However, the smallness of \(\epsilon\) can be exploited for a faster algorithm, which does not eliminate exponential growth in \(T\), but we are at least able to reduce bounding constants with lowered computational cost.

## 3 Small-\(\epsilon\) Asymptotic Analysis & Implementation

Let's consider the parameterization with \(\epsilon\) being small enough that effects of order \(\epsilon^{2}\) can be dropped or grouped in with round-off error. In this setting, we can construct a solution in a power series form and then truncate terms order \(\epsilon^{2}\) and higher. In other words, our expansion has a base term equal to the explicitly computable myopic solution (obtained for the case \(\delta=0\)), and an order-\(\epsilon\) correction term. Computation of the order-\(\epsilon\) correction is more involved, but computational costs and runtime are minimal.

### Expansion of \(\lambda_{t}\)

We write the formal expression for \(\lambda_{t}\),

\[\lambda_{t}=\epsilon\widetilde{\lambda}_{t}\;,\]

which we insert into (6) and (7) to obtain the following equations,

\[X_{t} =X_{t-1}-\frac{1}{q(\vec{S}_{t},P_{t})}\Psi_{t}^{-1}\widetilde{ \lambda}_{t} \tag{16}\] \[\widetilde{\lambda}_{t} =\delta\mathbb{E}_{t}\widetilde{\lambda}_{t+1}-\frac{\gamma}{ \epsilon}P_{t}\left(\frac{1}{\gamma}P_{t}^{-1}\Psi_{t}\mu-X_{t}\right)\] (17) \[=\delta\mathbb{E}_{t}\widetilde{\lambda}_{t+1}-\frac{\gamma}{ \epsilon}P_{t}\left(\frac{1}{\gamma}P_{t}^{-1}\Psi_{t}\mu-X_{t-1}\right)- \frac{\gamma}{\epsilon q(\vec{S}_{t},P_{t})}P_{t}\Psi_{t}^{-1}\widetilde{ \lambda}_{t}\;.\]
Rearranging (17) we have the following stabilized equation,

\[\widetilde{\lambda}_{t}=\left(\epsilon I+\frac{\gamma}{q(\vec{S}_{t},P_{t})}P_{t} \Psi_{t}^{-1}\right)^{-1}\left(\delta\epsilon\mathbb{E}_{t}\widetilde{\lambda} _{t+1}-\gamma P_{t}\left(\frac{1}{\gamma}P_{t}^{-1}\Psi_{t}\mu-X_{t-1}\right) \right)\, \tag{18}\]

for which it is straightforward to check that the convergence proof of Theorem 1 still applies. We write the following formal expansion,

\[\widetilde{\lambda}_{t}=\widetilde{\lambda}_{t}^{[0]}+\epsilon\widetilde{ \lambda}_{t}^{[1]}+\epsilon^{2}\widetilde{\lambda}_{t}^{[2]}+\dots\,\]

which we insert into (18) to obtain the following recursive expressions for the expansion's terms,

\[\widetilde{\lambda}_{t}^{[0]} =-\left(\epsilon I+\frac{\gamma}{q(\vec{S}_{t},P_{t})}P_{t}\Psi_{ t}^{-1}\right)^{-1}\left(\Psi_{t}\mu-\gamma P_{t}X_{t-1}\right) \tag{19}\] \[\widetilde{\lambda}_{t}^{[t]} =\delta\left(\epsilon I+\frac{\gamma}{q(\vec{S}_{t},P_{t})}P_{t} \Psi_{t}^{-1}\right)^{-1}\mathbb{E}_{t}\widetilde{\lambda}_{t+1}^{[i-1]}\]

for \(i=1,2,3,\dots\). This expansion is in powers of \(\epsilon\) and can be truncated for a good approximation of the solution to (16).

**Remark 2**.: _When we approximate \(\widetilde{\lambda}_{t}\) with lower-order terms \(\widetilde{\lambda}_{t}^{[0]}\) and \(\widetilde{\lambda}_{t}^{[1]}\), we can simplify the expressions in (19) so that they do not depend on \(\epsilon\)_

\[\widetilde{\lambda}_{t}^{[0]} =-\left(\frac{\gamma}{q(\vec{S}_{t},P_{t})}P_{t}\Psi_{t}^{-1} \right)^{-1}\left(\Psi_{t}\mu-\gamma P_{t}X_{t-1}\right) \tag{20}\] \[\widetilde{\lambda}_{t}^{[1]} =\delta\left(\frac{\gamma}{q(\vec{S}_{t},P_{t})}P_{t}\Psi_{t}^{- 1}\right)^{-1}\mathbb{E}_{t}\widetilde{\lambda}_{t+1}^{[0]}+\left(\frac{ \gamma}{q(\vec{S}_{t},P_{t})}P_{t}\Psi_{t}^{-1}\right)^{-2}\left(\Psi_{t}\mu- \gamma P_{t}X_{t-1}\right)\,\]

_so that \(\widetilde{\lambda}_{t}=\widetilde{\lambda}_{t}^{[0]}+\epsilon\widetilde{ \lambda}_{t}^{[1]}+\mathcal{O}(\epsilon^{2})\).The expansion in (20) is different from (19) because it has reduced the base term to the naive policy that goes straight to the aim portfolio, thus leaving it to the correction terms to compensate for transaction costs._

Inserting the \(\widetilde{\lambda}_{t}\) expansion into (16), we suspect the following lower-order expansion for the state process has a Big-Oh error as follows,

\[X_{t}=X_{t-1}-\frac{1}{q(\vec{S}_{t},P_{t})}\Psi_{t}^{-1}\left(\widetilde{ \lambda}_{t}^{[0]}+\epsilon\widetilde{\lambda}_{t}^{[1]}\right)+\mathcal{O}( \epsilon^{2}). \tag{21}\]

These order \(\mathcal{O}(\epsilon^{2})\) errors are in fact true, based on the following proposition:

**Proposition 2**.: _Assume Condition 1, Condition 2 and Condition 3. The order-\(\epsilon\) approximation \(\widetilde{\lambda}_{t}^{[0]}+\epsilon\widetilde{\lambda}_{t}^{[1]}\) has an error, which is_

\[\underset{t\leq T}{\sup}\mathbb{E}\left\|\widetilde{\lambda}_{t}^{[0]}+ \epsilon\widetilde{\lambda}_{t}^{[1]}-\widetilde{\lambda}_{t}^{\star}\right\| =\mathcal{O}\left(\epsilon\Delta(\epsilon)\exp\Big{(}\frac{\Omega T}{1- \Delta(\epsilon)}\Big{)}\mathbb{E}\underset{t\leq T}{\sup}\|\widetilde{ \lambda}_{t}^{[1]}\|\right)\,\]

_where \(\widetilde{\lambda}^{\star}\) denotes the solution to (16) and (17)._

Proof.: (see Appendix.)
```
Initialize: \(\theta^{(1)}\sim\mathcal{N}(0,0.01)\). for k = 1 to MAX_ITERdo  Initialize: \(P_{0}\), \(\vec{S}_{0}\) and \(X_{0}\).  Set: \(S_{0}=\text{diag}(\vec{S}_{0})\). for t = 1 to T do  ### Update the MGARCH state: \(Z_{t}=R_{t}-\mu\) \(\Sigma_{t}=CC^{\top}+A\Sigma_{t-1}A^{\top}+BZ_{t}Z_{t}^{\top}B^{\top}\) \(S_{t}=S_{t-1}(1+R_{t})\) \(\vec{S}_{t}=(S_{t}^{1},S_{t}^{2},...,S_{t}^{n})^{T}\) \(\Psi_{t}=\text{diag}(\vec{S}_{t})\) \(P_{t}=\Psi_{t}\Sigma_{t}\Psi_{t}\) ### Update the portfolio: \(\widetilde{\lambda}_{t}^{[0]}=-\left(\epsilon I+\frac{\gamma}{q(\vec{S}_{t}, P_{t})}P_{t}\Psi_{t}^{-1}\right)^{-1}\left(\Psi_{t}\mu-\gamma P_{t}X_{t-1}\right)\) \(\widetilde{\lambda}_{t}^{[1]}=\delta\left(\epsilon I+\frac{\gamma}{q(\vec{S}_ {t},P_{t})}P_{t}\Psi_{t}^{-1}\right)^{-1}\varphi(X_{t-1},\vec{S}_{t},P_{t}; \theta^{(k)})\) \(X_{t}=X_{t-1}-\frac{1}{q(\vec{S}_{t},P_{t})}\Psi_{t}^{-1}\Big{(}\widetilde{ \lambda}_{t}^{[0]}+\epsilon\widetilde{\lambda}_{t}^{[1]}\Big{)}\) endfor \(\text{loss}(\theta)=\frac{\sum_{t=1}^{T}\left\|\varphi(X_{t-1},\vec{S}_{t},P_ {t};\theta)-\widetilde{\lambda}_{t+1}^{[0]}\right\|^{2}}{T-1}\) \(\theta^{*}=\arg\min_{\theta}\text{loss}(\theta)\) \(\theta^{(k+1)}=\alpha\theta^{*}+(1-\alpha)\theta^{(k)}\) endfor
```

**Algorithm 1** Small-\(\epsilon\) Neural Network Fixed-Point Algorithm for MGARCH, with Learning Rate \(\alpha\in(0,1]\)

### Neural Network Algorithm

Let function \(\varphi:\mathbb{R}^{d}\times\mathbb{R}^{d}\times\mathbb{R}^{d\times d}\to \mathbb{R}^{n}\) be from a class of NN functionals with sigmoidal activation. The term \(\widetilde{\lambda}_{t}^{[1]}\) is approximated as \(\delta\left(\epsilon I+\frac{\gamma}{q(\vec{S}_{t},P_{t})}P_{t}\Psi_{t}^{-1} \right)^{-1}\varphi(X_{t-1},\vec{S}_{t},P_{t};\theta)\), where the optimal parameter \(\theta\) is found from

\[\min_{\theta}\sum_{t=1}^{T}\mathbb{E}\left\|\varphi(X_{t-1},\vec{S}_{t},P_{t}; \theta)-\mathbb{E}_{t}[\widetilde{\lambda}_{t+1}^{[0]}]\mathbb{I}_{t<T} \right\|^{2}\.\]

Note that now we are considering a NN architecture that remains constant through time, which is a considerable simplification from the NN architecture used in the proof of Theorem 2. The tradeoff in making this simplification is faster computation time. Algorithm 1 gives the implementation of the scheme in (13) using the lower-order expansion in (21) with this NN approximation of \(\delta\mathbb{E}_{t}\widetilde{\lambda}_{t+1}^{[0]}\). Fig. 2 shows the flow of Alg. 1 at each moment \(t\). We empirically observed that after training, our NN policy can work in real-time (specifically, each NN inference takes around 1.5 ms on a commodity laptop).

In our analysis of Algorithm 1's portfolio, we will compare to the purely myopic strategy,

\[X_{t}^{\text{myopic}}=X_{t-1}^{\text{myopic}}+\frac{1}{\epsilon q(\vec{S}_{ t},P_{t})}\Psi_{t}^{-1}\left(I+\widetilde{P}_{t}\Psi_{t}^{-1}\right)^{-1} \left(\mu-\gamma P_{t}X_{t-1}^{\text{myopic}}\right). \tag{22}\]

The myopic strategy shares the same objective as the RL strategy but for \(\delta=0\). Note that because our objective function contains a quadratic transaction cost term, the typical buy-and-hold portfolios will create a large transaction cost at the beginning, leading to a negative total wealth return. Therefore, we do not compare our method with the typical buy-and-hold portfolios.
## 4 Sector ETFs: an 11-Dimensional Example

We verify the small-\(\epsilon\) asymptotic analysis by allocating 11 sector ETFs, including iShares U.S. Real Estate ETF (IYR), Materials Select Sector SPDR Fund (XLB), Energy Select Sector SPDR Fund (XLE), Financial Select Sector SPDR Fund (XLF), Industrial Select Sector SPDR Fund (XLI), Technology Select Sector SPDR Fund (XLK), Consumer Staples Select Sector SPDR Fund (XLP), Utilities Select Sector SPDR Fund (XLU), Health Care Select Sector SPDR Fund (XLV), Consumer Discretionary Select Sector SPDR Fund (XLY), and Vanguard Communication Services ETF (VOX). The chosen ETFs are a good representation of a cross-section of the U.S. stock market returns. However, our method is applicable to any set of stocks. Training of Algorithm 1 was implemented on synthetic data sampled from an MGARCH model that was estimated on historical data of the 11 ETFs. We will show the implementation results of Algorithm 1 on both simulated testing data, and on historical out-of-sample data that was not used to estimate the MGARCH model. Overall, our algorithm practically performs well. Additionally, it is not necessary that the rigorous conditions of Sec. 2 have to hold, as we do not enforce Condition 2's boundedness of \(h\).

### Setup

#### 4.1.1 Dataset

We downloaded the adjusted closing price of the 11-sector ETFs in the year 2010 to the year 2019 from Yahoo Finance ([44]). The data was split into 10 folds for training, with each fold having 5 years of historical data. The starting and ending dates for the folds are shown in Table 1.

#### 4.1.2 Estimating Parameters

Given a sequence of historical price \(\{S_{t}^{i}\}\) for the \(i^{th}\) ETF, the return rate \(R_{t}^{i}\) at each time \(t\) can be found by:

\[R_{t}^{i}=\frac{S_{t+1}^{i}-S_{t}^{i}}{S_{t}^{i}}.\]

Figure 2: Flow of Alg. 1
Defining \(R_{t}=[R_{t}^{1},R_{t}^{2},...,R_{t}^{11}]^{\top}\) and \(\bar{R}\) the mean of \(R_{t}\) over \(t\), the initial covariance matrix \(\Sigma_{0}\) is calculated as:

\[\Sigma_{0}=\frac{1}{L}\sum_{t=1}^{L}(R_{t}-\bar{R})(R_{t}-\bar{R})^{\top},\]

where \(L\) is the length of the sequential historical data.

\(A\), \(B\), and \(C\) in (2) were estimated using Broyden Fletcher Goldfarb Shanno (BFGS) algorithm ([45]). Specifically, we first burned in the initial covariance matrix \(\Sigma_{0}\) into (2). We then used \(\sum_{t\geq 0}||\Sigma_{t+1}-\Sigma_{t}||^{2}\) with matrix norm as the loss function and BFGS as the minimizer to find \(A\), \(B\), \(C\), and \(\Sigma\) to minimize the loss function. The following equilibrium state is achieved when the loss function is the smallest:

\[\Sigma=CC^{\top}+A\Sigma A^{\top}+B\Sigma B^{\top}. \tag{23}\]

For the estimation of \(\mu\), we used the eigen-portfolio approach in [46]. The parameters (i.e., \(A\), \(B\), \(C\), \(\Sigma_{0}\), and \(\mu\)) were estimated for each fold. Fig. 3 shows an example of the historical ETF prices and the simulated ETF prices using the estimated parameters.

#### 4.1.3 Architecture of Neural Network

The utilized neural network (NN) is composed of fully connected layers. The input contains the portfolio \(X_{t-1}\) that has 11 elements, plus the covariance matrix of the dollar-returns \(P_{t}=\Psi_{t}\Sigma_{t}\Psi_{t}\) whose dimension is \(11\times 11\) and also the expected value of returns \(\mu^{\top}\Psi_{t}\) which is also 11-dimensional. Therefore, the total dimension is 143. The hidden layer size was determined by considering both the training time and the NN performance. When using more complex NN architecture, we observed that there was no obvious improvement in performance while the training time increased considerably. On the other hand, when using even simpler NN architectures, we observed that the deep RL algorithm suffered from under-fitting problems. Therefore, we utilized four hidden layers, each of which contains 400 neurons. The output of the NN corresponds to \(\varphi(\cdot,\cdot,\cdot;\theta)\) in Algorithm 1, which is 11-dimensional. The activation function is Tanh. The architecture of the utilized neural network (NN) is shown in Table 5 in the appendix shows the detail of the architecture. The programming language is Python 3. Tensorflow and Keras were utilized as the deep-learning library. Keras is built on top of Tensorflow.

### Simulation on Synthetic Data

At moment \(t+1\), a noise vector \(Z_{t+1}\) was generated from a Gaussian distribution \(\mathcal{N}(0,\Sigma_{t})\), where \(\Sigma_{t}\) is the covariance matrix at the previous moment, the return rate \(R_{t+1}\) was determined by (1), and the covariance matrix \(\Sigma_{t+1}\) was found by (2). For each fold, we generated one sequence of synthetic data to train the NN and the other 200 sequences to test the performance of the NN. We show how the total wealth grows with time for RL and myopic strategies. We also show how the fund manager should invest the money into the stock market with time using RL and myopic strategies, respectively. Specifically, wealth \(W_{t}\) at time \(t\) was calculated by

\[W_{t}=W_{t-1}+X_{t-1}^{\top}(\vec{S}_{t}-\vec{S}_{t-1})-\frac{\epsilon q(\vec {S}_{t},P_{t})}{2}a_{t}^{\top}\Psi_{t}a_{t}\,\]

\begin{table}
\begin{tabular}{c c c|c c c} \hline \hline
**Fold** & **From** & **To** & **Fold** & **From** & **To** \\ \hline
1 & Jan. 2010 & Oct. 2014 & 2 & Jul. 2010 & Apr. 2015 \\
3 & Jan. 2011 & Oct. 2015 & 4 & Jul. 2011 & Apr. 2016 \\
5 & Feb. 2012 & Nov. 2016 & 6 & Aug. 2012 & May 2017 \\
7 & Mar. 2013 & Nov. 2017 & 8 & Sep. 2013 & Jun. 2018 \\
9 & Feb. 2014 & Dec. 2018 & 10 & Aug. 2014 & Jun. 2019 \\ \hline \hline \end{tabular}
\end{table}
Table 1: Details of Each Fold
with \(W_{0}=100\). We ran the trained NN and the myopic strategy (as given by (22)) on the test data for each fold and calculated the average. The training epoch was 100 for each fold. The optimizer was Adam. The loss function was mean-squared-error loss (MSE). In this case, we took \(\epsilon=0.003\), for which we were able empirically to verify Condition 4. The results are shown in Fig. 4.

Table 2 shows the value added of the two methods in each fold on the synthetic data. The RL strategy of Algorithm 1 consistently outperforms the myopic strategy given by (22) by showing an annual increase of 1.8% in total value. The results are scalable to any \(W_{0}\). Therefore, our RL approach will show significant outperformance in the absolute value of total assets added when the given capital (i.e., \(W_{0}\)) is large. Note in the table that we are emphasizing the $-value of returns, as this is a better measure of fund performance as per the reasoning of [47].

The amount of money invested in the stock market is calculated by

\[I_{t}=I_{t-1}+a_{t-1}^{\top}\vec{S}_{t-1}\;,\]

where \(I_{0}=0\) because all money is in cash at the beginning. Note that since the ETF prices vary with time, it is likely that at the end of the day, the total money invested in the stock market may exceed \(W_{0}\). Fig. 5 shows the evolution of \(I_{t}\) for both RL and myopic strategies. The two strategies show different investment speeds. The RL is investing faster than the myopic strategy. Fig. 6 shows the evolution of the RL and myopic portfolio allocations with time. In the beginning, the portfolio allocations are at 0 for both RL and myopic. As time passes, the portfolios increase, but the slope becomes smaller. Eventually, the portfolios will be attracted to a stationary state. It is worth mentioning that the investment portfolio charts are not for comparison purposes. Instead, they are used to illustrate how our RL approach and the myopic approach are applied and to ensure that the two approaches have reasonable behaviors.

Figure 3: X-axis shows the day. Y-axis shows the ETF prices.
### Simulation on Historical Data

For historical out-of-sample data, \(Z_{t}\) was found by

\[Z_{t}=R_{t}-\mu\]

Figure 4: Average wealth (200 trajectories) of the RL strategy (solid) and the myopic strategy (dot). The X-axis shows the day. Y-axis shows the total wealth.
where \(R_{t}\) is the return rate of the real historical data at each time. For each fold, the NN was trained with the data of that fold and tested out-of-sample on the historical data of the following six months. The training epoch was 100. The optimizer was Adam. And the loss function was MSE. In this case, we took \(\epsilon=0.01\), and again we were able empirically to verify Condition 4. The result is shown in Fig. 7. For most situations, the RL outperforms the myopic strategy, but not as significantly as it did on synthetic data. This is because the real historical market data contains considerable uncertainty that might cause the diminished performance of the RL strategy in testing. Table 3 shows annualized total asset value added by using the two approaches on historical market data. Among the cases, the RL strategy shows more added value than the myopic approach. The RL method is able to consistently outperform the myopic by about 11bps. This over-performance is not tremendous but certainly makes evident that RL can improve out-of-sample performance for this problem. The results are scalable: if the fund manager has 100 million dollars, using the RL strategy can help her gain an extra $120,000 compared to the gain realized from the myopic strategy.

Fig. 8 shows how the money is invested with time. Fig. 9 shows how the portfolios change. It bears repeating that the practical purpose of this paper's optimization is to optimally move a large amount of new capital into the market, i.e., bring this new capital into the fund. By showing the cumulative amount of money invested we can see how much of the investment goal has been accomplished. Figures 5 and 8 illustrate how the two strategies behave as they are moving the new capital into stocks. Overall, the portfolios show similar behavior to what we observed in the synthetic data case, namely, that the RL strategy seeks to move the new capital into the ETFs faster than the myopic.
for the weaker significance when testing out-of-sample is because of the following. First, Fig. 4 shows persistent over-performance by RL because it is an average of 200 trajectories. It is possible that if we had abundantly more out-of-sample data, we could see averages playing out, in which case we might see a stronger out-of-sample performance by RL. Second, there may be an out-of-sample model risk, i.e., that the model changes or is misspecified in the out-of-sample test. In this case, the training is aimed at learning an optimum that becomes sub-optimal when applied in the out-of-sample test.

Figure 5: Money invested in the stock market at a different time for the RL strategy (solid) and the myopic strategy (dot). The X-axis shows the day. Y-axis shows the money invested in the stock market.
Ultimately, we cannot say much more than this, but the different folds of data that we have shown do give us some sense for out-of-sample variation for both RL and the myopic portfolios.

Figure 6: Portfolios of the RL strategy (solid) and the myopic strategy (dot) at different times. The X-axis shows the day. Y-axis shows the portfolio.
Figure 7: Total wealth of the RL strategy (solid) and the myopic strategy (dot) at different time. X-axis shows the date (month - day). Y-axis shows the wealth.
Figure 8: Money invested in the stock market at different times for the RL strategy (solid) and the myopic strategy (dot). The X-axis shows the date (month - day). Y-axis shows the money invested in the stock market.
### The Objective Value \(V\) on Historical Data

The main goal of this paper is to maximize the objective value function (4), we, therefore, ran the RL and myopic algorithms on the 10 folds historical market data and calculated the relevant terms in (4). Specifically, we show the

Figure 9: Portfolio value of the RL strategy (solid) and the myopic strategy (dot) with time. The X-axis shows the date (month - day). Y-axis shows the portfolio.
average results of \(V\), total transaction cost \(\sum_{t=0}^{T}\delta^{t}\frac{\epsilon q(\vec{S}_{t},P_{t})}{2}a_{t}^{\top}\Psi_{ t}a_{t}\), and risk penalty \(\sum_{t=0}^{T}\delta^{t}\frac{\epsilon}{2}X_{t}^{\top}P_{t}X_{t}\) with different \(\epsilon\). We observed that when \(\epsilon\) is relatively small (e.g., \(0.00001\)), our RL algorithm considers an aggressive investment to maximize (4) (i.e., increasing the inventory \(X_{t}\) quicker and thus taking a relatively higher total transaction cost and risk penalty). When \(\epsilon\) is relatively large (e.g., \(0.01\)), our RL algorithm adopts a conservative investment to maximize (4) (i.e., increasing \(X_{t}\) slower and thus taking a relatively lower total transaction cost and risk penalty). Fig. 10 shows how the inventory of IYR ETF changes for different \(\epsilon\) by using the RL and myopic strategies. The inventories of other ETFs show similar behavior. From the figure, when \(\epsilon=0.01\), our RL tends to slowly increase \(X_{t}\). When \(\epsilon=0.00001\), our RL tends to increase \(X_{t}\) quickly. The explanation for such behavior is that when \(\epsilon\) is small (meaning that the transaction cost can be negligible), \(f(a_{t},X_{t},\vec{S}_{t},P_{t})\) in (4) can be approximated as:

\[f(a_{t},X_{t},\vec{S}_{t},P_{t})\approx\underbrace{\mu^{\top}\Psi_{t}X_{t}}_{ \text{investment return}}-\underbrace{\frac{\gamma}{2}X_{t}^{\top}P_{t}X_{t}}_{ \text{risk penalty}} \tag{24}\]

which has an equilibrium point at \(X_{t}=\frac{1}{\gamma}P_{t}^{-1}\mu^{\top}\Psi_{t}\) for every \(t\). Therefore, to maximize (4) with the approximation (24), our RL algorithm needs to increase \(X_{t}\) to the equilibrium point within a short time and thus performs an aggressive investment. If \(\epsilon\) is large (meaning that the transaction cost cannot be negligible), to avoid high transaction costs, our RL can only slowly increase \(X_{t}\). We show the average results of final objective-function value \(V\), total transaction costs, and total risk penalty over the 10 folds historical market data in Table 4 for different \(\epsilon\). The total transaction cost and risk penalty for \(\epsilon=0.00001\) are higher than the total transaction cost and risk penalty for \(\epsilon=0.01\), which confirms our explanation.

Our RL algorithm outperforms the myopic strategy for both aggressive and conservative investment cases by showing a higher \(V\) according to Table 4. For the aggressive investment case (\(\epsilon=0.00001\)), our RL algorithm behaves less aggressively than the myopic strategy by having a lower total transaction cost and risk penalty. On the other hand, for the conservative investment case (\(\epsilon=0.01\)), our RL algorithm behaves more aggressively than the myopic strategy by taking a higher total transaction cost and risk penalty. However, for both cases, our RL algorithm shows a higher value for the objective \(V\) than the myopic strategy, as shown in Table 4.

Fig. 11 shows the average performance of our RL and the myopic strategies with time on the 10-fold historical market datasets. The average wealth returns using RL and myopic strategies are close. However, our RL approach attains a higher (i.e., better) value for the objective than the myopic strategy, meaning that our RL achieves a better trade-off
between investment return, transaction cost, and risk penalty. Therefore, our RL outperforms the myopic strategy.

### Discussion of Our RL Approach

From the experimental results, we have observed that our RL approach has a higher objective value V given in (4) (i.e., the main goal of this paper) than the myopic strategy on the historical data. In other words, our RL approach balances the wealth return, risk penalty, and transaction cost better than the myopic by showing a 3.85% better for \(\epsilon=0.01\) and 11.64% better for \(\epsilon=0.00001\). On synthetic data, the RL method returns an extra 1-2% additional annual percentage return over the myopic (out of a total of roughly 10% annual return). This demonstrates that our method outperforms the myopic, as this extra 1-2% annually is significant when measuring the growth of assets. For real data, we are able to show consistent over-performance of the RL method by about 11bps. This over-performance certainly makes it evident RL's improved out-of-sample performance. Our RL approach also shows a faster investment speed with working in real time. Therefore, our approach is a valid effective method.

## 5 Conclusion

We have implemented a reinforcement learning algorithm to solve the mean-variance preferences of (4) with an MGARCH model and transaction costs of order \(\epsilon\), where \(\epsilon>0\) is a small parameter. Our method addresses issues of algorithm convergence and computation time by using an \(\epsilon\) expansion to guide the network to the solution. While similar methods use double (dueling) networks to stabilize convergence, we only need to use one network with our expansion. Our method is stable and can work in real time. The resulting portfolios show good performance in simulated tests.

## Appendix A Proofs

**Proof of Proposition 1**: We apply the Sherman-Woodbury-Morrison formula,

\[\left(I+\widetilde{P}_{t}\Psi_{t}^{-1}\right)^{-1}=\Psi_{t}\widetilde{P}_{t}^ {-1}-\Psi_{t}\widetilde{P}_{t}^{-1}\left(I+\Psi_{t}\widetilde{P}_{t}^{-1} \right)^{-1}\Psi_{t}\widetilde{P}_{t}^{-1}. \tag{25}\]

The first term in (25) is bounded as follows,

\[\|\Psi_{t}\widetilde{P}_{t}^{-1}\|=\frac{\epsilon}{\gamma}q(\widetilde{S}_{t}, P_{t})\|P_{t}^{-1}\Psi_{t}\|\leq\frac{\epsilon}{\gamma}\chi\|h\|_{\infty}\|P_{t}^{- 1}\|\leq\frac{\epsilon\chi\|h\|_{\infty}}{\gamma\underline{\varsigma}^{2}c}\,\]

\begin{table}
\begin{tabular}{|c|c c c c|} \hline \multirow{2}{*}{\(\epsilon\)} & \multicolumn{4}{c|}{\(V\)} \\ \(\epsilon\) & RL & Myopic & Difference & Difference / RL \\ \hline \(0.01\) & 0.0234 & 0.0225 & 0.0009 & 3.85\% \\ \(0.00001\) & 0.0885 & 0.0782 & 0.0103 & 11.64\% \\ \hline \multirow{4}{*}{\(0.01\)} & \multicolumn{4}{c|}{Total Transaction Cost} \\  & RL & Myopic & Difference & Difference / RL \\ \hline \(0.01\) & 0.0013 & 0.0011 & 0.0002 & 15.38\% \\ \(0.00001\) & 0.2562 & 0.2600 & -0.0038 & -1.50\% \\ \hline \multirow{4}{*}{\(0.01\)} & \multicolumn{4}{c|}{Total Risk Penalty} \\  & RL & Myopic & Difference & Difference / RL \\ \hline \(0.01\) & 0.0012 & 0.0011 & 0.0001 & 8.33\% \\ \(0.00001\) & 0.2399 & 0.2436 & -0.0037 & -1.54\% \\ \hline \end{tabular}
\end{table}
Table 4: Total Value \(V\), Transaction Costs, and Risk Penalties for Different \(\epsilon\).
Figure 11: Results on real data. Top: \(\epsilon=0.00001\). Bottom: \(\epsilon=0.01\).
where we have bounded \(\|P_{t}^{-1}\|\) in the following way,

\[\|P_{t}^{-1}\|\leq\|\Psi_{t}^{-1}\|^{2}\cdot\|\Sigma_{t}^{-1}\|\leq \frac{1}{\underline{\mathrm{s}}^{2}\inf_{\|v\|=1}v^{\top}CC^{\top}v}=\frac{1}{ \underline{\mathrm{s}}^{2}c}\,\]

where \(c>0\) is the lower bound defined in Condition 1. The second term in (25) is bounded as follows,

\[\left\|\Psi_{t}\widetilde{P}_{t}^{-1}\left(I+\Psi_{t}\widetilde{P}_{t}^{-1} \right)^{-1}\Psi_{t}\widetilde{P}_{t}^{-1}\right\|\leq\left\|\Psi_{t} \widetilde{P}_{t}^{-1}\right\|^{2}\left\|\left(I+\Psi_{t}\widetilde{P}_{t}^{-1 }\right)^{-1}\right\|\leq\left(\frac{\epsilon X\|h\|_{\infty}}{\gamma \underline{\mathrm{s}}^{2}c}\right)^{2}\left\|\left(I+\Psi_{t}\widetilde{P}_{ t}^{-1}\right)^{-1}\right\|\.\]

Now, taking the norm of (25), applying a triangle inequality to the right-hand side, and inserting these bounds, we have

\[\left\|\left(I+\Psi_{t}\widetilde{P}_{t}^{-1}\right)^{-1}\right\|\leq\frac{ \epsilon X\|h\|_{\infty}}{\gamma\underline{\mathrm{s}}^{2}c}+\left(\frac{ \epsilon X\|h\|_{\infty}}{\gamma\underline{\mathrm{s}}^{2}c}\right)^{2}\left\| \left(I+\Psi_{t}\widetilde{P}_{t}^{-1}\right)^{-1}\right\|\.\]

If \(\epsilon X\|h\|_{\infty}<\gamma\underline{\mathrm{s}}^{2}c\), then we can re-arrange to obtain the bound in (9).

**Proof of Theorem 1**: The following lemma is useful to have when proving convergence in Theorem 1.

**Lemma 1**.: _Let \(A\) and \(B\) be two symmetric positive definite matrices. Then \(\|(A+B)^{-1}\|\leq\min(\|A^{-1}\|,\|B^{-1}\|)\)._

Proof.: By positive definiteness of each matrix, we have

\[\inf_{\|v\|=1}v^{\top}(A+B)v\geq\inf_{\|v\|=1}v^{\top}Av+\inf_{\|v\|=1}v^{\top }Bv=\frac{1}{\|A^{-1}\|}+\frac{1}{\|B^{-1}\|}\geq\frac{1}{\min(\|A^{-1}\|,\|B^ {-1}\|)}\.\]

Thus, we have

\[\|(A+B)^{-1}\|=\frac{1}{\inf_{\|v\|=1}v^{\top}(A+B)v}\leq\min(\|A^{-1}\|,\|B^ {-1}\|)\,\]

which proves the second statement. 

From inspection of (11), the usefulness of the following lemma should be evident:

**Lemma 2**.: _If we assume Condition 1, Condition 2 and Condition 3, then in (11) there are the following bounds in coefficients,_

\[\left\|\frac{1}{\epsilon q(\widetilde{S_{t}},P_{t})}\Psi_{t}^{-1}\right\| \leq\frac{1}{\epsilon\underline{\mathrm{s}}} \tag{26}\] \[\left\|\left(I+\widetilde{P}_{t}\Psi_{t}^{-1}\right)^{-1}P_{t}\right\| \leq\frac{\epsilon}{\gamma}\chi\|h\|_{\infty}\,\]

_where \(\underline{\mathrm{s}}\) and \(\|h\|_{\infty}\) are constants defined in Condition 2 and \(X\) is the bound defined in Condition 3._

Proof.:
1. The first bound in (26) is a clear consequence of Condition 2 and Condition 3.
2. To prove the second bound in (26), we start as follows, \[\left(I+\widetilde{P}_{t}\Psi_{t}^{-1}\right)^{-1}P_{t}=\frac{\epsilon}{\gamma }q(\widetilde{S_{t}},P_{t})\left(\frac{\epsilon}{\gamma}q(\widetilde{S_{t}},P_ {t})P_{t}^{-1}+\Psi_{t}^{-1}\right)^{-1}\,\] where \(P_{t}\) is invertible because we have assumed Condition 1 putting a lower bound on the eigenvalues. However, we don't have a lower bound on \(\|P_{t}^{-1}\|\) because the eigenvalues of \(P_{t}\) do not have a fixed upper bound, but from Condition 2 we have \(\|\Psi_{t}\|\leq\|h\|_{\infty}\), and thus we use Lemma 1 to get \[\frac{\epsilon}{\gamma}q(\widetilde{S_{t}},P_{t})\left\|\left(\frac{\epsilon}{ \gamma}q(\widetilde{S_{t}},P_{t})P_{t}^{-1}+\Psi_{t}^{-1}\right)^{-1}\right\| \leq\frac{\epsilon}{\gamma}q(\widetilde{S_{t}},P_{t})\|h\|_{\infty}\,\] which proves the second bound in (26).
We are now ready to prove Theorem 1. From (11) and using \(\Omega\) defined in (12), we obtain the following inequality,

\[\|\lambda_{t}^{(k+1)}-\lambda_{t}^{(k)}\|\leq\delta\|(I+\widetilde{P} _{t}\Psi_{t}^{-1})^{-1}\|\cdot\mathbb{E}_{t}\|\lambda_{t+1}^{(k+1)}-\lambda_{t+ 1}^{(k)}\|\] \[+\|\gamma(I+\widetilde{P}_{t}\Psi_{t}^{-1})^{-1}P_{t}\|\sum_{t^{ \prime}=1}^{t-1}\frac{\|\Psi_{t^{\prime}}^{-1}\|}{\epsilon q(\widetilde{S}_{t^ {\prime}},P_{t^{\prime}})}\|\lambda_{t^{\prime}}^{(k)}-\lambda_{t^{\prime}}^{( k-1)}\|\] \[\leq\Delta(\epsilon)\mathbb{E}_{t}\|\lambda_{t+1}^{(k+1)}-\lambda _{t+1}^{(k)}\|+\Omega\sum_{t^{\prime}=1}^{t-1}\|\lambda_{t^{\prime}}^{(k)}- \lambda_{t^{\prime}}^{(k-1)}\|\, \tag{27}\]

where from Condition 4 we have \(\left\|\delta\left(I+\widetilde{P}_{t}\Psi_{t}^{-1}\right)^{-1}\right\|\leq \Delta(\epsilon)<1\) and from Lemma 2 it follows that

\[\sup_{t,t^{\prime}}\frac{\|\Psi_{t^{\prime}}^{-1}\|}{\epsilon q(\widetilde{S} _{t^{\prime}},P_{t^{\prime}})}\left\|\left(I+\widetilde{P}_{t}\Psi_{t}^{-1} \right)^{-1}P_{t}\right\|\leq\Omega\.\]

If we consider a finite-time version of the problem where \(\lambda_{T+1}^{(k)}\equiv 0\) for all \(k\) (i.e., no control on \(X_{t}^{(k)}=X_{T}^{(k)}\) for \(t>T\)), then we can prove convergence of the iteration scheme. Denote the expected iteration error as

\[\mathcal{E}_{t}^{(k)}=\mathbb{E}\|\lambda_{t}^{(k)}-\lambda_{t}^{(k-1)}\|\.\]

Denoting the vector of errors as

\[\mathcal{E}_{1:T}^{(k)}=\left(\mathcal{E}_{T}^{(k)},\mathcal{E}_{T-1}^{(k)}, \ldots,\mathcal{E}_{1}^{(k)}\right)^{\top},\]

the inequality in (27) can be expressed with the following matrix/vector system,

\[\mathcal{E}_{1:T}^{(k+1)}\leq H_{T}\mathcal{E}_{1:T}^{(k+1)}+N_{T}\mathcal{E} _{1:T}^{(k)} \tag{28}\]

where \(H_{T}\) the \(T\times T\) nilpotent matrix whose entry at row \(i\) column \(j\) is

\[H_{T}^{ij}=\Delta(\epsilon)\times\begin{cases}1\text{ if }i=j+1\\ 0\text{ otherwise,}\end{cases}\]

which has norm \(\|H_{T}^{k}\|=\Delta^{k}(\epsilon)\) for \(k<T\) and \(\|H_{T}^{k}\|=0\) for \(k>T\), and where \(N_{T}\) is the \(T\times T\) nilpotent matrix whose entry at row \(i\) column \(j\) is

\[N_{T}^{ij}=\Omega\times\begin{cases}1\text{ if }i<j\\ 0\text{ otherwise,}\end{cases}\]

for which \(N_{T}^{k}=0\) for \(k\geq T\). Therefore, by combining the error bound in (28) we have

\[\|\mathcal{E}_{1:T}^{(k+1)}\|\leq\left\|H_{T}\mathcal{E}_{1:T}^{(k+1)}\right\| +\left\|N_{T}\mathcal{E}_{1:T}^{(k)}\right\|\leq\Delta(\epsilon)\left\| \mathcal{E}_{1:T}^{(k+1)}\right\|+\left\|N_{T}\mathcal{E}_{1:T}^{(k)}\right\|\,\]

and after rearranging and taking \(k\)-many iterations, we have

\[\|\mathcal{E}_{1:T}^{(k+1)}\|\leq\left(\frac{1}{1-\Delta(\epsilon)}\right)^{k }\|N_{T}^{k}\mathcal{E}_{1:T}^{(1)}\|\leq\left(\frac{1}{1-\Delta(\epsilon)} \right)^{k}\|N_{T}^{k}\|\|\mathcal{E}_{1:T}^{(1)}\|. \tag{29}\]

Now we take a moment to derive a bound on the normal of \(N_{T}\): comparing sums with integrals we see the following,

\[\|N_{T}^{k}e_{T}\|_{1}\leq\Omega^{k}\sum_{t_{1}=1}^{T-1}\sum_{t_{2}=1}^{t_{1}- 1}\cdots\sum_{t_{k}=1}^{t_{k-1}-1}1\leq\Omega^{k}\int_{0}^{T}dt_{1}\int_{0}^{ t_{1}}dt_{2}\cdots\int_{0}^{t_{k-1}}dt_{k}=\frac{\Omega^{k}T^{k}}{k!}\,\]
for \(k<T\) where \(e_{T}\in\mathbb{R}^{T}\) is the \(T^{th}\) canonical basis vector, and this bound gives us the general bound

\[\|N_{T}^{k}\|\leq\sqrt{T}\|N_{T}^{k}e_{T}\|_{1}\leq\frac{\Omega^{k}T^{k+1/2}}{k! }\mathbb{1}_{k<T}\qquad\forall k. \tag{30}\]

We apply the bound in (30) to the right-hand side in (29), and we find that \(\|\mathcal{E}_{1:T}^{(k+1)}\|=0\) for all \(k\geq T\) thus proving convergence to a unique fixed point, hence proving the statement of Theorem 1.

**Proof of Theorem 2**: Given the compact set \(\mathbf{K}\) from (14), by the universal approximation theorem ([25]), for each iteration in (13) there is a NN for which parameters can be chosen such that

\[\sup_{t\leq T}\mathbb{E}\Big{[}\,\Big{\|}\lambda(t,X_{t-1}^{(k)},\vec{S}_{t},P _{t};\theta^{(k)})-Y_{t}^{(k-1)}(\theta^{(k)})\Big{\|}\,\mathbb{1}_{\{(X_{t-1} ^{(k)},\vec{S}_{t},P_{t})\in\mathbf{K}\}}\Big{]}\leq\varepsilon^{(k)}\,\]

where \(\varepsilon^{(k)}\) is the arbitrarily small error that can be decreased by increasing the NN's hyperparameters. In the same manner as the process given by (8), we define the iterated process to be

\[\lambda_{t}^{(k)}=\lambda(t,X_{t-1}^{(k)},\vec{S}_{t},P_{t},\theta^{(k)})\.\]

Denote the error as

\[\mathcal{E}_{t}^{(k)}=\mathbb{E}\,\Big{\|}\lambda_{t}^{(k)}-\lambda_{t}^{*} \Big{\|}\,\]

where \(\lambda_{t}^{*}\) is the limiting fixed-point of (8). Using the pair of iterative equations in (8), we can include the NN error to obtain the following recursive bound,

\[\mathcal{E}_{t}^{(k)}\leq\frac{\mathbb{E}\,\Big{[}\|\gamma(I+\widehat{P}_{t} \Psi_{t}^{-1})^{-1}P_{t}\|\cdot\|X_{t-1}^{(k)}-X_{t-1}^{*}\|\Big{]}+(\varepsilon ^{(k)}+2\eta)}{1-\Delta(\epsilon)}\.\]

Then, proceeding similarly to the proof of Theorem 1 in Appendix, we denote the vector of all errors as \(\mathcal{E}_{1:T}^{(k)}=\left(\mathcal{E}_{T}^{(k)},\mathcal{E}_{T-1}^{(k)}, \ldots,\mathcal{E}_{1}^{(k)}\right)^{\top}\), and like we had in (28), a bound on these errors can be expressed with the following matrix/vector system,

\[\mathcal{E}_{1:T}^{(k+1)} \leq\frac{1}{1-\Delta(\epsilon)}\left(N_{T}\mathcal{E}_{1:T}^{(k) }+\sup_{\ell}(\varepsilon^{(\ell)}+2\eta)\mathbf{1}\right)\] \[\leq\left(\frac{1}{1-\Delta(\epsilon)}\right)^{k}N_{T}^{k} \mathcal{E}_{1:T}^{(1)}+\sup_{\ell}(\varepsilon^{(\ell)}+2\eta)\sum_{i=0}^{k- 1}\left(\frac{1}{1-\Delta(\epsilon)}\right)^{i+1}N_{T}^{i}\mathbf{1}\,\]

where \(H_{T}\) and \(N_{T}\) are the same matrices that were defined in Appendix, and where \(\mathbf{1}\) denotes the vector in \(\mathbb{R}^{T}\) of all 1's. Thus,

\[\|\mathcal{E}_{1:T}^{(k+1)}\|\leq\left(\frac{1}{1-\Delta(\epsilon )}\right)^{k}\|N_{T}^{k}\mathcal{E}_{1:T}^{(1)}\|+\sup_{\ell}(\varepsilon^{( \ell)}+2\eta)\sum_{i=0}^{k-1}\left(\frac{1}{1-\Delta(\epsilon)}\right)^{i+1} \|N_{T}^{i}\mathbf{1}\|\] \[\leq\left(\frac{1}{1-\Delta(\epsilon)}\right)^{k}\|N_{T}^{k}\| \|\mathcal{E}_{1:T}^{(1)}\|+\sqrt{T}\left(\sup_{\ell}\varepsilon^{(\ell)}+2 \eta\right)\sum_{i=0}^{k-1}\left(\frac{1}{1-\Delta(\epsilon)}\right)^{i+1}\|N _{T}^{i}\|\] \[=\left(\frac{1}{1-\Delta(\epsilon)}\right)^{k}\frac{\Omega^{k}T^{k +1/2}}{k!}\,\mathbb{1}_{k<T}\|\mathcal{E}_{1:T}^{(1)}\|+\sqrt{T}\left(\sup_{ \ell}\varepsilon^{(\ell)}+2\eta\right)\sum_{i=0}^{k-1}\left(\frac{1}{1-\Delta (\epsilon)}\right)^{i+1}\frac{\Omega^{i}T^{i+1/2}}{i!}\,\mathbb{1}_{i<T}\,\]

for large \(k\geq T\) we have

\[\|\mathcal{E}_{1:T}^{(k+1)}\|=\mathcal{O}\left(\left(\sup_{\ell}\varepsilon^{( \ell)}+2\eta\right)\exp\left(\frac{\Omega T}{1-\Delta(\epsilon)}\right)\right)\,\]

which is the statement of Theorem 2.
**Proof of Proposition 2**: Denote the state process from the order-\(\epsilon\) approximation as

\[X_{t}^{[1,2]}=X_{t-1}^{[1,2]}-\sum_{t^{\prime}=1}^{t}\frac{1}{q(\widetilde{S_{t^ {\prime}}}_{t^{\prime}}P_{t^{\prime}})}\Psi_{t^{\prime}}^{-1}(\tilde{\lambda}_{t ^{\prime}}^{[0]}+\epsilon\widetilde{\lambda}_{t}^{[1]})\,\]

for \(t=1,2,3,\ldots,T\). Inserting the expansion into (17) results in the order-\(\epsilon\) approximation error

\[\widetilde{\lambda}_{t}^{[0]}+\epsilon\widetilde{\lambda}_{t}^{[ 1]}-\left(\epsilon I+\frac{\gamma}{q(\widetilde{S_{t}},P_{t})}P_{t}\Psi_{t}^{- 1}\right)^{-1}\left(\delta\epsilon\mathbb{E}_{t}\left[\widetilde{\lambda}_{t+1 }^{[0]}+\epsilon\widetilde{\lambda}_{t+1}^{[1]}\right]-\gamma P_{t}\left( \frac{1}{\gamma}P_{t}^{-1}\Psi_{t}\mu-X_{t-1}^{[1,2]}\right)\right)\] \[=-\left(\epsilon I+\frac{\gamma}{q(\widetilde{S_{t}},P_{t})}P_{t} \Psi_{t}^{-1}\right)^{-1}\left(\delta\epsilon^{2}\mathbb{E}_{t}\widetilde{ \lambda}_{t+1}^{[1]}\right)\,\]

where the equality uses the expressions of for \(\widetilde{\lambda}_{t}^{[0]}\) and \(\widetilde{\lambda}_{t}^{[1]}\) given in (19). Using Lemma 2, we have

\[\left\|\left(\epsilon I+\frac{\gamma}{q(\widetilde{S_{t}},P_{t})}P_{t}\Psi_{t} ^{-1}\right)^{-1}\right\|=\frac{\Delta(\epsilon)}{\epsilon}\.\]

Therefore the differential error given above is of order \(\mathcal{O}\left(\epsilon\Delta(\epsilon)\mathbb{E}\sup_{t\leq T}\|\widetilde{ \lambda}_{t}^{[1]}\|\right)\). Thus, similar to Theorem 2, we have a big-oh bound on error,

\[\sup_{t\leq T}\mathbb{E}\left\|\widetilde{\lambda}_{t}^{[0]}+\epsilon \widetilde{\lambda}_{t}^{[1]}-\widetilde{\lambda}_{t}^{\star}\right\|= \mathcal{O}\left(\epsilon\Delta(\epsilon)\exp\Big{(}\frac{\Omega T}{1- \Delta(\epsilon)}\Big{)}\mathbb{E}\sup_{t\leq T}\|\widetilde{\lambda}_{t}^{[1 ]}\|\right)\,\]

where \(\widetilde{\lambda}^{\star}\) denotes the solution from (16) and (17).
