Receding Horizon Optimization of Large Trade Orders

Received April 9, 2021, accepted April 20, 2021, date of publication April 26, 2021, date of current version May 3, 2021. Digital Object Identifier 10.1109/ACCESS.2021.3075700 Receding Horizon Optimization of Large Trade Orders GERARDO ALFONSO 1,2, A. DANIEL CARNERERO 2, DANIEL R. RAMIREZ 2, AND TEODORO ALAMO 2 1Shenwan and Hongyuan Securities Co., Ltd. (SWS), Shanghai 200031, China 2Departamento de Ingeniería de Sistemas y Automática, Universidad de Sevilla, 41092 Sevilla, Spain Corresponding author: Daniel R. Ramirez ([email protected]) This work was supported by the Agencia Estatal de Investigación (AEI), Spain, under Grant PID2019-106212RB-C41/AEI/10.13039/501100011033. ABSTRACT The optimal execution of stock trades is a relevant and interesting problem as it is key in maximizing profits and reducing risks when investing in the stock market. In the case of large orders, the problem becomes even more complex as the impact of the order in the market has to be taken into account. The usual solution is to split large orders into a set of smaller suborders that must be executed within a prescribed time window. This leads to the problem of deciding when in the time window execute each suborder. There are popular ways of executing the trading of these split orders like those which try to track the ‘‘Time Weighted Average Price’’ and the ‘‘Volume Weighted Average Price’’, usually called TWAP and VWAP orders. This paper presents a strategy to optimize the splitting of large trade orders over a given time window. The strategy is based on the solution of an optimization problem that is applied following a receding horizon approach. This approach reduces the impact of prediction errors due to the uncertain market dynamics, by using new values of the price time series as they are available as time goes on. Suborder size constraints are taken into account in both market and limit orders. The strategy relies on price and traded volume forecast but it is independent of the prediction technique used. The performance index weighs not only the financial cost of the suborders, but also the impact on the market and the forecasting accuracy. A tailored optimization algorithm is proposed for efficiently solving the corresponding optimization problem. Most of the computations of the algorithm can be parallelized. Finally, the proposed approach has been tested through a case study composed by stocks of the Chinese A-share market. INDEX TERMS Algorithmic trading, receding horizon optimization, large stock orders, limit orders, TWAP, VWAP. I. INTRODUCTION Stock trading is becoming an increasingly complex field as investors try to maximize their profits and reduce their risks with increased emphasis in recent years on optimal execution. After an investor has decided to buy or sell stocks, there are many options on how to execute this trade. In large stock orders, it is very frequent to split orders into smaller orders in an attempt not to impact the market. For instance, a buy order to purchase a large amount of stock could push the price of the stock up if executed in one single block. The market impact of that order can be potentially reduced by splitting that single order into smaller orders and executing them over The associate editor coordinating the review of this manuscript and approving it for publication was Hao Wang . time rather than in one go. There are many ways in which this can be accomplished. One of the simplest approaches is called ‘‘Time Weighted Average Price’’, commonly referred to as TWAP. A TWAP order splits the order in blocks of shares of the same size that are then executed at regular time intervals. A more sophisticated trade is a ‘‘Volume Weighted Average Price’’, commonly referred to as VWAP. This type of order is very common as it represents around 50% of all the institutional investors’ trading [1]. In a VWAP order, rather than slicing the original order in smaller trades of equal size (equal number of shares), a forecasted traded volume for the desired interval is estimated. The size of each (sliced) transaction is proportional to the forecasted volume for the corresponding time bucket. Thus, a critical step in this approach is to be able to generate volume forecasts. Besides these two popular VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ 63865 G. Alfonso et al.: Receding Horizon Optimization of Large Trade Orders ways of trading sliced orders, there are basically two types of trading orders: market orders and limit orders. A market order is an order to be executed immediately at the prevailing market price. The focus on this type of order is speed rather than price optimization. Limit orders, on the other hand, focus on price efficiency. A limit order necessarily has an associated price over (below) which the buy (sell) order cannot be executed. There is a tradeoff between speed of execution and price optimization. Moreover, in the case of limit orders, there is no certainty of execution as it is possible to not have enough time buckets in the execution period for which the price meets the limit. Thus, in these orders, success is measured not only by the price attained, but also by the percentage of the order that has been executed. There is relatively limited existing literature on the topic of optimal execution of split trade orders using learning techniques, with more papers covering stock forecasting by means of different techniques like neural networks [2] or deep learning [3], support vector machines [4], [5], adaptive line combiners [6] or local data-based techniques [7]. There are however some interesting articles in the field of trade execution optimization such as for instance [8]. In this article, the authors proposed a genetic algorithm to optimize a limit order book used for price formation in an artificial stock market. Genetic algorithms are also used in [9] to generate trading strategies, not based on forecasting, that are back-tested against historical data of the Australian Stock Exchange. Recurrent neural networks have been used in [10] to predict price-flip events in limit order books by classifying sequences of observations of the book depths and market orders. Forecasting the traded volume has also been used in [11] to improve the execution of VWAP orders. That is, by forecasting the traded volume one can track the VWAP price matching it at the end of the chosen time window. On the other hand, [12] derives analytical solutions of a static optimal execution strategy of a VWAP trade, in which the optimal execution strategy can be calculated by an iteration of a single variable optimization, rather than by a multivariable optimization. In that work, the market is modelled using non-anticipating and Brownian motion processes. Particle swarm optimization has been used by [13] in a high-frequency trading system based on moving averages, used to determine the trading sequence that maximizes the net returns over a series of consecutive time steps. This optimization technique has also been used by [14] to train a kernel-based nonlinear predictor that was also applied to forecast the VWAP price in the Shanghai market. Optimal control methods have been also considered for generating the sequence of suborders in split large orders. For instance, [15] presented an optimal VWAP algorithm based on the linear quadratic regulator (LQR) subject to limits in the size of the suborders. The linear model for the stock prices is based on Brownian motion, a type of model that has also been used by [16] to find the stock prices also in the context of VWAP operations. Hamilton- Jacobi-Belmann methods and in general variational calculus have been used in [17]–[19] based on Brownian motion and random walks models. A recent paper [20] has studied optimal VWAP strategies using unconstrained optimization on models based on the assumption that the stock prices can be modelled as martingales and the traded volume as autoregressive processes. Time series methods have been used in [21], where an autoregressive fractionally integrated moving average model is used to forecast intraday trading volumes in the Chinese equity market, and its application to VWAP tracking, obtaining better results than static approaches. VWAP tracking has also been tackled in [22], where an interesting combination of historical averages and SVM to forecast intraday trading volumes in the gold and S&P 500 futures markets was used. Another interesting approach to forecasting trading volumes in a VWAP tracking context has been presented in [23]. In this paper, the authors used the fast Fourier transform algorithm to identify the periodic and the non-periodic part of the trading volume, using historical values of 50 stocks contained in the Shanghai 50 stock index. A similar approach is followed in [24]. In this paper, we propose a data-based method that relies on price and volume forecasts, and the use of dynamic optimization over a finite horizon, to obtain optimal sequences of suborders to fulfil large trade orders. The optimal sequence is computed by minimizing a cost index, in which several terms are taking into account. Besides the price, the trading impact factor is considered along with a term related to forecasting accuracy. The technique is based on solving an optimization problem each time bucket of the time window in which the order has to be executed. Thus, at each time bucket, a complete sequence of suborders for the remaining time window is obtained, but only the first component of the sequence is effectively used. This strategy is similar to the feedback receding horizon or predictive control strategy used in automatic control [25], but with a shrinking prediction horizon. Using such an approach allows us to split the order into different suborders each of them computed with the most recent available information (and thus with improved forecastings). Therefore, this approach reduces the effects of the prediction errors that are bound to arise, especially if the prediction horizon is large. The optimization is carried out by means of a tailored, efficient algorithm that can solve the optimization problem very fast. Furthermore, most of its computations can be done in parallel if necessary. This computational efficiency can be very useful if high-frequency trading is considered, or if learning or simulation-based strategies are used to tune the hyperparameters of the algorithm. Regarding the type of orders, both market and limit orders are considered in the paper and treated in a unified way. Finally, as a case study, stocks from the onshore Chinese A-share market are used to show the performance of the proposed strategy. The rest of the paper is outlined as follows: Section II describes the problem statement for each type of order. Section III presents the optimization algorithm. Section IV 63866 VOLUME 9, 2021 G. Alfonso et al.: Receding Horizon Optimization of Large Trade Orders presents the results of applying the proposed technique to the case study. Finally, section Vpresents the conclusions. II. PROBLEM STATEMENT The objective of this paper is to design a strategy to execute large stock orders that, in order to limit their impact on the market, have to be split into a number of smaller suborders. Thus, an order to buy M∈ <+shares will be executed by splitting the order into up to Nsuborders m(t+k)∈ <+, with k∈[1,N], such that the sum of all the suborders is equal to Mor, in the case of limit orders, the closest possible to M. Note that in practice, the number of shares to be traded must be a natural number, so the number of shares in the suborders that are sent to the market are rounded to the nearest smaller integer, that is bm(t+k)c. Given that the number of shares considered in these orders is quite large, often hundreds of thousands or millions of shares per order, the difference in financial cost between considering real or integer-valued suborders is negligible. However, from a computational point of view, the difference is quite significant, being much more efficient an algorithm based on real-valued suborders. The proposed strategycomputes the splitting of the original order in an optimal way, that is, achieving the most convenient price. The strategy is based on forecasting both the price of the stock and its total traded volume over the time window defined by N, i.e., from t+1 to t+N. The price and trade volume forecasts are used to compute a performance index that should be optimized. Thus the strategy will rely on an optimization problem in which a performance index, denoted as VN, is optimized. The formulation of VNmust weigh certain aspects of concern, such as the total cost of the order, but also the impact on the market of each of the suborders. Besides that, given that forecasts are used to compute the performance index, a term penalizing the potential degradation of forecasts with the prediction horizon is also included in VN. Let ˆp(t+k|t) be the price forecast for t+k, ˆv(t+k|t) the total traded volume forecast1for t+kand m(t+k|t) the number of shares to be bought at t+k, meaning the notation t+k|tthat these values are computed at time t. Then the proposed performance index is: VN(mN(t),ˆ pN(t),ˆ vN(t)) = N X k=1aˆp(t+k|t) +αm(t+k|t) ˆv(t+k|t)β!m(t+k|t) +µ N X k=1 σkm(t+k|t),(1) where a=1 in buy orders or a= −1 in sell orders, the tuning parameters α,µ,σ,βare nonnegative with β≥1, and mN(t)=[m(t+1|t),...,m(t+N|t)], 1It is assumed that price and volume forecast are greater than zero. While this is met almost always for the price, the trade volume can be indeed zero, thus zero volume forecasts should be changed for an arbitrarily small value. ˆ pN(t)=ˆp(t+1|t),...,ˆp(t+N|t), ˆ vN(t)=ˆv(t+1|t),...,ˆv(t+N|t), the sequences of suborders, price and traded volume forecastings respectively. Note that the first term represents the forecasted amount or economic figure of the buy or sell operation, and that in this term an impact factor correction has been included. Impact factor correction represents the influence of the suborder volume in the price for that time bucket. Impact factors can be modelled as linear terms, like in [15] or be more elaborate like the exponential form used here based on [26]. The second term assigns a greater cost to suborders that are far in the future, because the prediction error grows with the prediction horizon. Note that for this effect the exponential weight σkmust have σ > 1. The proposed strategy aims to find the optimal sequence of suborders m∗ N(t)=[m∗(t+1|t),m∗(t+2|t),...,m∗(t+N|t)],(2) that minimizes the performance index VNover the time window N. However, depending on the type of order, some constraints have to be taken into account. In the case of market orders, the suborders must meet a certain size limit so that the impact in the market is limited, i.e., 0≤m(t+k|t)≤m,∀k∈[1,N], being m=0.1Ma typical value (i.e., a suborder cannot exceed a 10% of the whole order). On the other hand, in limit orders, the size limit will be the same, but only if the price forecast meets the price limit pl. Notice that, in the case of a buy limit order, ˆp(t+k|t)>plimplies that m(t+k|t)=0, that is, the buy suborder cannot be sent to the market. Similarly, in sell limit orders ˆp(t+k|t)<plimplies that m(t+k|t)=0. These situations can be addressed with the following constraint 0≤m(t+k|t)≤mk,∀k∈[1,N] (3) where mk=(0 if limit order & a(ˆp(t+k|t)−pl)>0 motherwise (4) On the other hand, the sum of all suborders should be, in principle, equal to the total number of shares to be traded (M) thus the equality constraint N X k=1 m(t+k|t)=M should be taken into account. In the case of limit orders, the possibility of not having enough time buckets in which the price limit is met, must be taken into consideration. To illustrate this consider that if every suborder must be at most 10% of the original order and the price forecast is only under the limit in 7 time buckets, the original order can only be completed at most at a 70%, which would be equal to VOLUME 9, 2021 63867 G. Alfonso et al.: Receding Horizon Optimization of Large Trade Orders the sum of the corresponding mkfor all the time window. That is, when the order cannot be fully executed, the largest percentage of fulfilment of the order is attained forcing N X k=1 m(t+k|t)= N X k=1 mk. Note that if N X k=1 mk≥M then the order can be fully executed. Thus, to be able to address both situations, we include the following constraint that maximizes the degree of fulfilment of the order, N X k=1 m(t+k|t)=Mc,(5) with Mc=min(M, N X k=1 mk).(6) Note that this constraint will work for either market or limit orders. Once all the necessary constraints have been formulated, the optimal sequence of suborders will be obtained by solving the optimization problem: m∗ N(t)=arg min mN(t)VN(mN(t),ˆ pN(t),ˆ vN(t)) s.t. (3) and (5).(7) After problem (7) is solved, one could apply the entire optimal sequence m∗(t) executing the suborders m∗(t+ 1|t),...,m∗(t+N|t) at the corresponding time buckets. This approach suffers from two related issues: first, as the index kgrows, the forecastings ˆp(t+k|t) and ˆv(t+k|t) have a higher prediction error. Thus the suborders m(t+k) will rely on progressively more inaccurate forecastings and thus the computed optimal value will differ from the ideal optimal value that could be computed if the real values of p(t+k) and v(t+k) were known in advance. Second, as time advances, new real values of p(t+k) and v(t+k) are available, but they are not used. These new values could be used to obtain better predictions of the remaining time window, allowing to compute optimal suborders that would be closer to their ideal values, i.e., those obtained without prediction errors. Thus, in this paper it is proposed to use a feedback receding horizon optimization strategy, typical of predictive control techniques [25], in which at each time bucket tproblem (7) is solved to obtain m∗ N(t) sending only m∗(t+1|t) to the market and discarding the rest of the sequence m∗ N(t). To apply this receding horizon strategy it is necessary to take into account that the time window length shrinks with time, thus an initial time window length, denoted as Np, must be considered. In this way, initially N=Np, and then it will be decreasing at each time step. Furthermore, constraint (5) must reflect the fact that Mwould be decreasing as suborders are being sent to the market. Thus, an initial size order MIwill be considered, i.e., for the first time M=MI, and then updated as necessary. Besides the aforementioned changes, the receding horizon strategy can help in mitigating the negative effect of price prediction error in limit orders. This error can affect in two ways, one is that the price limit is met by the predicted price, whereas the real price fails to meet the limit. The other is the reverse, that is, the predicted price does not meet the limit, but the real one does. The solution to these issues is to compute the optimal sequence of suborders for the case that ˆp(t+1|t) fails to meet pland for the case that it meets pl, and then, at time t+1, apply the one that corresponds to the value of p(t+1). For the first case, it is not really necessary to perform any computation as it is obvious that m∗(t+1|t)=0. For the second case, problem (7) must be solved with the bound m1=mregardless the value of ˆp(t+1|t). So, this modification can be easily implemented by just changing (4) to mk=     0 if limit order & a(ˆp(t+k|t)−pl)>0 & k>1 motherwise (8) and sending to the market m∗(t+1|t) only if p(t+1) meets pl. Algorithm 1summarizes the proposed strategy. Algorithm 1 Optimal Execution of Large Orders Require: MI,Np,a, α, β, µ, σ, m. 1: M←MI. 2: N←Np. 3: repeat 4: Compute the forecasts ˆ pN(t) and ˆ vN(t). 5: Compute the bounds mkas in (8). 6: Compute Mcas in (6). 7: Solve (7) to obtain m∗ N(t). 8: Wait for next time bucket, that is, t+1. 9: Compute m∗ t+1=(0 if limit order & a(p(t+1)−pl)>0 m∗(t+1|t) otherwise 10: if m∗ t+1>0then 11: Send m∗ t+1to the market. 12: M←M−m∗ t+1. 13: end if 14: N←N−1. 15: until N=0 or M=0. Remark 1: The initial prediction horizon Npcan be quite long; thus the prediction errors can affect the performance of the proposed strategy. The receding strategy helps to mitigate such effects by using forecastings obtained with the most up to date information. Note also that the time window shrinks at each step, reducing the prediction horizon and allowing us to work with progressively better forecasts. Furthermore, 63868 VOLUME 9, 2021 G. Alfonso et al.: Receding Horizon Optimization of Large Trade Orders the receding horizon only sends to the market the first suborder in m∗ N(t), that mostly depends on short term forecastings, which are the most accurate. Finally, the forecasting errors in long prediction horizons are also compensated by the second term of (1), assigning more weight to suborders that are placed far in the prediction horizon depending on the accuracy of the forecasting method. The performance of the strategy depends, not only on the price and traded volume forecastings, but also on the choice of the performance function (1) hyperparameters. Regarding the tuning of such hyperparameters, αand βshould be chosen by the user so that they reflect the impact of each suborder in the stock price. This implies a modelling task that is previous to the decision of investing in a particular stock and that it is beyond the scope of this paper (see [26]). On the other hand, σand µcan be chosen considering the length of the prediction horizon and the long term accuracy of the forecasting method used to predict the price time series. Alternatively, the algorithm hyperparameters can be tuned by performing a search of optimal values running simulations based on historical data, choosing the values that optimize a certain performance metric, e.g., the total financial cost of executing the whole order. Finally, note that at each time bucket, problem (7) must be solved. Because of the high frequency of nowadays trading systems, it is clear that an efficient solving method for (7) is needed. Furthermore, having an efficient optimization algorithm is important if the optimization of the hyperparameters is carried out, as this implies solving (7) many times. III. OPTIMIZATION ALGORITHM The nonlinear optimization problem (7) is relatively easy to solve, as the performance index and constraints are convex. Note that the equality constraint (5) is the cause of most of the computational burden of problem (7). This constraint is the one that is coupling the Ndecision variables of problem (7). Without this constraint, the Ndecision variables optimization problem could be solved as Nindependent problems of just one decision variable. This would be much more efficient from a computational point of view, as these Nproblems could even be solved in a parallel way. In order to decouple the decision variables, we resort to a dual formulation. Consider, for a given λ≥0, the problem: J∗(λ)=min mN(t)VN(mN(t),ˆ pN(t),ˆ vN(t)) +λ Mc− N X k=1 mk! s.t. 0 ≤mk≤mk,∀k∈[1,N],(9) where, for notational convenience, mkdenotes m(t+k|t). This optimization problem does not have in general the same solution as (7), but as will be shown in the following, for a certain value of λits solution would be the same and can be obtained with a closed formula. Taking into account the definition on VNin (1) and that this problem is free of the constraint (5), the solution of (9) can be computed by solving for every k∈[1,N] the following problem: m∗ k(λ)=arg min mk Jk(mk, λ) (10) s.t. 0 ≤mk≤mk(11) where m∗ k(λ) denotes the optimal value of mkfor a given λ and Jk(mk, λ)=aˆp(t+k|t)+αmk ˆv(t+k|t)β +µσk−λ mk which in turn can be rewritten as Jk(mk, λ)=(dk−λ)mk+ekmβ+1 k with dk=aˆp(t+k|t)+µσk and ek=α ˆv(t+k|t)β≥0.(12) Note that each m∗ k(λ) for k=1,...,Ncan be computed independently, thus the Nproblems (10) can be solved in parallel. To find m∗ k(λ) consider the partial derivative of Jk(mk, λ): ∂Jk(mk, λ) ∂mk =dk−λ+(β+1)ekmβ k. From (11) and (12) we have that mkand ekare non-negative. This implies that the partial derivative is strictly positive if dk−λ > 0. Thus, in this case, Jk(mk, λ) is a monotonically growing function of mkand the minimum is attained at mk= 0. That is, m∗ k(λ)=0 if λ≤λ− k=dk.(13) On the other hand, if dk−λ≤0, then the zero of the derivative, i.e. the minimizer of the unconstrained problem is attained at λ−dk (β+1)ek1 β .(14) We notice that, in this case, the minimizer of the unconstrained problem coincides with the minimizer of the constrained problem only if λ−dk (β+1)ek1 β ≤mk, This is equivalent to λ≤λ+ k=mβ k(β+1)ek+dk.(15) VOLUME 9, 2021 63869 G. Alfonso et al.: Receding Horizon Optimization of Large Trade Orders Thus, for all λ≥λ+ kthe constraint (11) determines the minimizer, i.e., m∗ k(λ)=mk. Summing up, the solution of (10) is m∗ k(λ)=         0 if λ≤λ− k λ−dk (β+1)ek1 βif λ∈[λ− k, λ+ k] mkif λ≥λ+ k, (16) where λ− kand λ+ kare given in (13) and (15). Note that the values of λproviding values different from 0 and mkare those in the interval [λmin, λmax], where λmin =min k∈[1,N]λ− k,(17) and λmax =max k∈[1,N]λ+ k.(18) Once the value of m∗ k(λ) has been obtained, consider M∗(λ) defined as M∗(λ)= N X k=1 m∗ k(λ).(19) In what follows, it will be shown that there exist λ∗∈ [λmin, λmax] such that M∗(λ∗)=Mc,(20) which implies that constraint (5) holds and the combined solutions m∗ k(λ∗) will be equal to that of (7), i.e., m∗ N(t)=[m∗ 1(λ∗),m∗ 2(λ∗),...,m∗ N(λ∗)].(21) The following property is key to show that λ∗exists. Property 1 (Extreme Values of m∗ kand M∗): The extreme values of m∗ k(λ) and M∗(λ) as functions of λare (m∗ k(λmin)=0,∀k∈[1,N] M∗(λmin)=0, and        m∗ k(λmax)=mk,∀k∈[1,N] M∗(λmax)= N X k=1 mk≥Mc. Proof: These values stem directly from (16) and the definitions of λmin and λmax. Taking into account that λmin ≤ λ− k, it is clear that m∗ k(λmin)=0, for all k. This in turn, makes that, by definition, M∗(λmin)=0. On the other hand, λmax ≥λ+ k, which implies that m∗ k(λmax)=mk, for all k. Furthermore, M∗(λmax)= PN k=1mk, which due to (6) is greater or equal to Mc. Note that, taking into account (16) it is obvious that m∗ k(λ) is continuous, as at λ=λ− k,m∗ k(λ− k)=0 from both sides and at λ=λ+ kits value is m∗ k(λ+ k)=mk. TABLE 1. Stock dataset. The following levels of market capitalization were follow for stock classification: small (<80 bn RMB), mid (from 80 to 300 bn) and large (>300 bn). Moreover, consider the derivative of m∗ k(λ), which is ∂m∗ k(λ) ∂λ = 1 β((β+1)ek)1 β!(−dk+λ)1 β−1 if λ∈(λ− k, λ+ k) and 0 otherwise. Taking into account that ek≥0 and that, by definition, −dk+λ− k=0, and −dk+λ+ k≥ 0, it follows that the derivative is nonnegative, thus m∗ k(λ) is a nondecreasing function. Given that m∗ k(λ) is a continuous nondecreasing function and considering (19), it follows that M∗(λ) is also a continuous nondecreasing function. Furthermore, Mcis nonnegative and by property 1, Mc∈"0, N X k=1 mk#=[M∗(λmin),M∗(λmax)]. This implies that a certain value λ∗∈[λmin, λmax] exists, such that M∗(λ∗)=Mc. The value λ∗that makes (20) and, therefore, (5) hold, is computed by performing a simple bisection search in the interval [λmin, λmax]. Once this λ∗is found, the optimal value of the suborders, that is m∗ k(λ∗), will already be computed, and the solution of problem (7) will be (21). Remark 2: The proposed strategy can be easily adapted to other forms of stock orders which require other constraints, like the so-called ‘‘max’’ orders in which the suborder size is also limited by a fraction of the forecasted traded volume for each time bucket. In these cases, the multiple inequality constraints should be reduced to one using a simple redundant constraint elimination. IV. CASE STUDY: CHINESE STOCK MARKET The proposed approach has been validated using stocks from the Chinese Stock Market. Every working day has two continuous trading sessions, each one two hours long. Thus, the objective of the case study will be to split large orders over a 120 minute session, being the time bucket 1 minute long. To validate the proposed strategy, nine stocks from the Chinese Stock Market have been chosen, splitting the choices between small, mid or large-capitalization stocks. Table 1 show the names and tickers of the nine stocks ordered according to their capitalization. The validation consists on the optimization of market orders of 2,000,000 shares over a 63870 VOLUME 9, 2021 G. Alfonso et al.: Receding Horizon Optimization of Large Trade Orders FIGURE 1. Validation data set for 002040. FIGURE 2. Validation data set for 688012. FIGURE 3. Validation data set for 000333. validation set of 20 two hour sessions randomly chosen from the range 1/4/2017 to 12/29/2020 (except the newer stock 688012 which is 1/2/2020 to 12/29/2020). For conciseness, only three of the validation datasets are shown in figures 1 to 3. It can be seen that the sessions in each set show bearish, bullish and sideways trends. The proposed trading strategy relies on forecasts of prices and traded volumes, but it is not linked to any particular forecasting technique. Here, neural networks and the local data-based technique proposed in [7] have been considered. In order to forecast prices using either technique, the market state, denoted as zp(t) has been described as composed by the last 100 prices and the last traded volume,2i.e., zp(t)=[p(t),p(t−1),...,p(t−99),v(t)]. 2This paper is not focused on the forecasting task but on the optimization of the splitting of the large order. Better forecastings could be obtained by carefully testing and choosing other market state vectors, e.g., using technical indicators, but this task is outside the scope of this paper. On the other hand, to forecast the traded volume, the market state, denoted by zv(t) is formed in this case by zv(t)=[v(t),v(t−1),...,v(t−99),p(t)]. The components of both state vectors are lightly smoothed with an exponential moving average of 5 minutes (EMA5), that is, pEMA(t)=2 6p(t)+1−2 6pEMA(t−1), but the target prices to be forecasted are unsmoothed. The number of market states in the local data set used by the technique of [7] has been chosen to be 350 and the weighting factor γ=0, whereas in the case of the neural network predictors, a Multilayer Perceptron with one hidden layer of 40 neurons has been used. Note that in the case of the local data technique of [7] it is not necessary to train the predictor, whereas in the case of the neural network, 120 networks have been trained using the Levenberg-Marquardt rule, one for each time bucket in the time window. It is noteworthy that the prediction horizon, determined by the time window and time bucket, is quite long, being Np=120, thus it is quite challenging in relation to the prediction errors, that are bound to arise. Taking into account that both price and volume have to be forecasting using one of these techniques, four possible prediction models have been considered, one using only neural networks (denoted as NN-NN), other relying only on the local data approach of [7] (LD-LD), and two combining these two techniques, i.e., LD-NN that predicts prices using the local data approach and volumes with neural networks and the reverse case, i.e., NN-LD. We first analyze market orders in which the objective is to split the large stock order, over each of the sessions of the validation set, optimizing the buying price. Thus the proposed strategy must solve (7) with an initial prediction horizon of N=Np=120 that will shrink down to N=1 as the receding horizon optimization is applied through the entire time window. Other parameters of VNused in this case study are, a=1 (buy orders), α=0.001, β=1.1, µ=0.6, σ= 1.05. The average prices obtained when using the proposed approach over the validation set of each stock with each of the four prediction models are shown in table 2. It can be seen that no model can be chosen as the best for all the stocks, even if they are close in their results (especially for models LD-LD and LD-NN). In fact, when considering the size of the orders, it is clear that even minute differences in prices count, thus, in general, for each stock it would be necessary to pick the best available model. To validate the results, two well-known benchmark prices will be considered for comparison purposes. The first one is the VWAP price assuming knowledge of the traded volume (which has to be forecasted in practice). This price is computed as VWAP(t) =PNp k=1v(t+k)p(t+k) PNp k=1v(t+k).(22) VOLUME 9, 2021 63871 G. Alfonso et al.: Receding Horizon Optimization of Large Trade Orders TABLE 2. Average prices in the validation set (RMB) using the proposed approach with each of the prediction models. The best value for each stock highlighted in bold. TABLE 3. Validation of the results by comparing to benchmark prices (RMB). The second baseline strategy is the TWAP price, which is the average price of the stock over the time window, that is, TWAP(t)= Np X k=1 p(t+k) Np .(23) Furthermore, the price attained with the proposed strategy if historical data, i.e., not forecasted, are used will be also considered. This price cannot be achieved in practice, but it is useful to find the maximum performance that could be reached in the absence of prediction errors, and therefore properly calibrate their impact on the performance of the strategy. Table 3shows the average prices over the validation set using both historical and forecasted prices and volumes, along with the VWAP and TWAP prices. It can be seen that the proposed strategy obtains a lower price than VWAP or TWAP in all the validation sets.3Clearly, the prediction error in the price and volume forecasting reduces the difference, but even with forecasted prices and volumes the proposed strategy obtains lower prices than the baseline approaches. Regarding the capitalization of the stocks, it appears that the proposed strategy gets worse results in the case of small-cap stocks, but this could be related to the differences in the price per share, getting better results in stocks with a higher price per share. On the other hand, it is quite remarkable that for some of the stocks, even the prices obtained with historical data are very close to VWAP or TWAP, e.g. stocks with ticker 300239 and 600000. In these cases, there is clearly a very small margin 3Note that, in the simulation, the impact of each suborder in the market has not been taken into account, as this impact is very difficult to model (the term in the performance function (1) acts as a de-tuning factor and does not constitute an accurate model of market impact). Thus, in practice, the results could be even better provided that the parameters αand βare properly chosen and adapted to the market conditions. TABLE 4. Savings on the validation set over VWAP and TWAP when using historical or forecasted data (RMB). of improvement, meaning that in these cases it proved very difficult to beat the market on the validation set. There are other stocks, like 688012, in which the hypothetical benefit that could be obtained is greatly wasted by the prediction errors. Regarding this, it is noteworthy that the prediction horizon is quite large in this case study, forcing the forecast up to 120 steps ahead. This is much longer than usual, given the fact that most of the stock forecasting applications focus on one step ahead predictions. On the other hand, higher prediction errors do not always imply a worse result as long as the price trend is correctly predicted. The reason for this is that in the proposed strategy price forecasting accuracy is not what really allows better results, but the ability to forecast the position of the lower prices in the time window. Finally, even if the price difference attained in practice seems small, given the high number of shares traded in this type of orders, the overall benefit can be significant. To illustrate this, consider table 4in which the savings of the proposed strategy over VWAP and TWAP are shown for the case of using the real or forecast prices and volumes. It can be seen that the savings are substantial in most cases and very substantial in some cases like the stocks with tickers 688012 and 600519. Furthermore, considering the savings that can be obtained over the many large orders that are traded over a year, it is evident that even a modest improvement on the price can justify the use of strategies like the proposed in this paper. The proposed approach has also been validated when using limit orders with the case study. Table 5shows the average prices and the percentage of order execution in the validation set for all the stocks. Different price limits have been set for each session in the validation set according to the price levels so that the price limit is not very stringent but not trivial to meet (which would result, in practice, in a de facto market order). Regarding the benchmark price comparison, it is noteworthy that in the case of limit orders the TWAP and VWAP prices cannot be considered, as they do not consider the limit in the price. Instead, TWAP limit and VWAP limit orders have been used as benchmarks. These orders work in the same way as conventional TWAP and VWAP orders except they are forced to meet the price limit. Thus, a TWAP or VWAP limit order cannot be started until the price limit is met, and they have to be stopped whenever the current price does not meet the limit. If the limit is met 63872 VOLUME 9, 2021 G. Alfonso et al.: Receding Horizon Optimization of Large Trade Orders TABLE 5. Average prices and percentage of execution in the case of limit orders. TABLE 6. Limit orders, stock 002040. again, a new TWAP or VWAP limit order is started for the remaining quantity to be executed and set to be executed in the remaining time window. As a result, TWAP and VWAP limit orders have no guarantee of being executed completely and their success can be judged by the percentage of the original order that it is executed. Thus, when comparing the proposed strategy to these benchmarks, the price should be taken into account only if the percentage of execution is the same. Note, however, that the attained price must always meet the price limit or the order will be considered failed (and the brokerage firm forced to pay the excess over the limit). The results in table 5show that the proposed strategy manages to reach a higher percentage of order execution than TWAP or VWAP limit orders while meeting the price limit. It is clear that the forecasting errors induce negative effects in both the price and execution percentage attained, but even using forecasted data the proposed strategy beats both TWAP and VWAP limit orders while meeting the price limit. Limit orders are more complex than market orders, thus it is interesting to examine in more detail the results shown in table 5. For conciseness, the results of only three stocks are shown in tables 6to 8(these are the same stocks whose dataset is shown figures 1to 3). It can be seen that the proposed approach usually achieves higher percentages of execution, although prediction errors can have a great impact on this percentage, e.g., session 18 in stock 002040 has 100% execution when using historical data but only 10% when using forecasted data. Even though, the latter result nearly doubles that of TWAP and VWAP limit orders. Another fact that can be observed is that there are sessions in which all the approaches are able to complete the orders. In these situations, the proposed approach using historical data achieves a better price, but when using forecasted data the price can be worse (e.g., session 7 in stock 002040 or session 6 in stock 688012 in which the results are better than the TWAP limit but worse than the VWAP limit). This is to be expected, as no approach can be the best in all situations in a task as difficult as trying to beat the market. Nevertheless, the proposed approach show on average better results (as illustrated in tables 5to 8). Furthermore, there are sessions in which the proposed approach, even with forecasted data, obtains a better price and a better percentage of execution (e.g., session 9 in stock 002040 or session 17 in stock 688012). On the other hand, the price limit can be so stringent that even with historical data is not possible to execute completely the order. This can be seen in sessions 4 and 19 of stock 002040 or in session 3 of stock 000333. Note that in those cases the proposed approach reached higher execution percentages than the benchmarks. Finally, regarding the computational burden of the optimization algorithm, the average time to solve the optimization problem in a serial Matlab implementation on a budget mobile processor Core i5-8265U clocked at 1.6 GHz is just 1.2 milliseconds. This is much lower than the time bucket but VOLUME 9, 2021 63873