A Monte Carlo method for optimal portfolio executions
aa r X i v : . [ q -f i n . T R ] D ec A Monte Carlo method for optimal portfolio executions
Nico Achtsis and Dirk NuyensJuly 24, 2018
Abstract
We treat the problem of mean-variance optimal execution in markets with limited liquidity andvarying volatility. When the market parameters are assumed constant, an analytical solution existsfor the optimal trading rate. In general however, this problem leads to a non-linear Hamilton–Jacobi–Bellman PDE, which has to be solved numerically. Since solving such a PDE is a complexprocedure, Almgren [2012] mentions a sub-optimal control that can be used as an approximation.This strategy assumes the market parameters are constant, and hence takes the analytical solutionfrom the stationary problem, but updates the strategy each time the market parameters change. Itis called the rolling horizon strategy (RHS), because it is essentially a continuously updated staticcontrol with contracting horizon. It is easy to extend to the multi-asset case as we will show. Inthis paper, we propose a rolling horizon Monte Carlo algorithm (RHMC). Our method chooses atrading rate based on simulations using a sub-optimal control. The potential upside of this methodis that our proposed RHMC method not only uses current market information, such as the RHS,but also uses simulations to infer future market behaviour as well. Our new method is naturallyformulated for the multi-asset case and allows the freedom to choose the structure of the stochasticdriver processes. The results indicate that our method can significantly outperform the RHS. Wealso provide some insights into the RHS, showing that it converges to the optimal solution for strongrisk-averse traders, at least in the setting of Almgren [2012].
Key words.
Quasi-Monte Carlo (QMC), Optimal asset execution, Optimal control problems.
AMS subject classifications.
This paper revolves around a fundamental part of algorithmic trading, namely, trade scheduling. Whenfacing the execution of a large block of assets over a fixed time interval, it is usually beneficial to split upthe order in several smaller blocks over the time interval to reduce market impact. Finding the optimumschedule requires balancing between market risk and liquidity risk. The former entails the risk of adverseprice moves of the assets that are traded, for this reason slower trading will lead to higher market risk.Liquidity risk on the other hand corresponds to the difference in the pretrade price before the orderis executed and the actual execution price. This difference is called slippage, and the accumulation ofthese gaps in prices will be larger when trading is faster. Consequently trade scheduling poses a dilemmabetween trading fast to eliminate market risk against trading slow to minimize slippage. We formulatethis so-called optimal execution problem in terms of a standard mean-variance optimization.The first market impact models were based on the discrete-time models constructed by Bertsimas and Lo[1998] and Almgren and Chriss [2001] and their continuous-time variants proposed by Almgren [2003].These models all separated the impact into two components: an instantaneous one affecting only theindividual trades that also triggered it, and a permanent effect that has an impact on all future trades.Research in market microstructure however suggests that market impact decays over time, as one cansee in the overview of, e.g., Eisler et al. [2012]. The first models to pick up on this fact are those ofObizhaeva and Wang [2013] and Potters and Bouchaud [2003]. The latter is an example of a limit orderbook model, meaning the authors model the dynamics of supply and demand in the order book to findthe optimal execution algorithm. This model was further developed by Alfonsi et al. [2010, 2012], inparticular to include nonlinear price impacts.An important aspect of modelling the market impact is consistency. Gatheral [2010] provides a goodoverview of what properties are desirable for such a model; e.g., one should not allow an algorithm to1anipulate the market through trading in order to profit. Huberman and Stanzl [2004] were among thefirst to point out that it is not sufficient to require the absence of arbitrage strategies in the usual sense.They illustrated that the feedback of trading strategies can lead to so-called price manipulation strategiesthat, when suitably rescaled and repeated, can create a weak form of arbitrage. Again, Gatheral [2010]explores the relation between impact functions and such weak forms of arbitrage. Furthermore, it wasshown by Alfonsi et al. [2012] that transaction-triggered price manipulation is possible in models that donot allow for price manipulation in the sense of Huberman and Stanzl [2004]. The papers by Alfonsi et al.[2010, 2012] provide models without such phenomena.The goal of this paper is to complement the work on single asset optimal execution schemes byRobert Almgren, in particular that of Almgren [2012]. Therein the author assumes a trader perceivesthe asset price plus an instantaneous effect based on the current trading speed. The problem of optimalexecution is there solved under the assumption that volatility and liquidity vary perfectly inversely(termed coordinated variation). In this paper we extend the model to multiple assets and do not needthe assumption of coordinated variation. As a minor side step, in Almgren [2012, Section 1.4] a simplifiedsolution called the rolling horizon strategy (RHS) is also proposed. This strategy entails that a stationarysolution is sought, which is updated during the trading period. In general this will not be optimal, butit is argued to be easy to implement while providing a reasonable solution. In this paper, we proposea rolling horizon Monte Carlo algorithm (RHMC) in which the trading rates are calculated based onsimulations using a sub-optimal control. In this way our new RHMC method not only uses currentmarket information, such as the RHS, but also uses simulations to infer future market behaviour.Almgren [2012] deals with optimal execution within a mean-variance framework, meaning that notonly the expected cost of trading is considered, but also the risk profile of the trader. Other papersusually deal with execution assuming a risk-neutral trader (and hence only look at the expected costincurred from trading). There are some papers which have already made multi-asset extensions to theAlmgren framework, see for instance Konishi [2002] and Sch¨oneborn [2011]. These papers however onlydeal with constant liquidity and volatility, and assume that there is no cross-liquidity effect betweenassets.In Section 2 we give an overview of the optimal execution problem. Under the assumption that marketparameters are constant, we derive the optimal strategy in Section 3. In Section 4 we look at strategieswhich allow dynamic market parameters. Section 4.1 contains two results on the RHS, showing whenthis strategy becomes optimal as well as its behaviour for risk-neutral traders. The main contribution isin Section 4.2, where we explain our rolling horizon Monte Carlo method (RHMC). This method is thentested numerically in Section 5. Finally we conclude in Section 6. h sec:EP i The execution problem consists in liquidating or acquiring n asset positions over a finite trading period[0 , T ], T < ∞ . We use the following notation: let x i ( t ) denote the number of i th shares that still needto be bought or sold at time t . Then x i (0) equals the number of shares that need to be traded at theinception of the program, where x i (0) < x i (0) > T with x i ( T ) = 0. We will denote the initialposition by the vector x . The trading speed, or first derivative of x i , will be denoted by v i . We canchoose either x i or v i as the control variable, since both determine the other. However, for numericalstability we will formulate our algorithms in terms of x . We use the notation x ′′ i to denote the secondderivative of x i with respect to time t . h sec:Model i We consider a probability space (Ω , F , P ) endowed with a filtration F = ( F ( t )) t ∈ [0 ,T ] which representsthe information structure available to the agent. We assume that F (0) is trivial and that F ( T ) = F .We also suppose that F satisfies the usual conditions of right-continuity and completeness (see, e.g.,Karatzas and Shreve [1991]). All components of the model will be defined on the filtered probabilityspace (Ω , F , F , P ). 2 isky asset The price of the risky assets follow an arithmetic Brownian motion S k ( t ) = S k (0) + Z t σ k ( s ) dW k ( s ) (1) h eq:S˙dynamics i with σ k ( t ) a (stochastic) function of time and dW k ( t ) dW ℓ ( t ) = ρ kℓ dt . The process W k ( t ) is a standard F ( t )-adapted Brownian motion. This model is known as Bachelier’s model and is widely used in theoptimal execution literature. See, among others, Alfonsi et al. [2012], Almgren [2012], Bertsimas and Lo[1998]. We prefer this model because of its simplicity, and standardised nature. We will write S ( t ) forthe vector containing S ( t ) to S n ( t ). Under the Bachelier model the dynamics for S ( t ) imply S ( T ) ∼ N S (0) , Z T Σ( s ) ds ! , where the matrix Σ( s ) consists of the elements σ k ( s ) σ ℓ ( s ) ρ kℓ . The natural assumption is made that thematrix Σ( s ) is positive definite. The Bachelier model could lead to negative asset values; however, since T is typically small, the probability of negative asset values is negligible. Price impact function
The price impact function gives the price change relative to the risky asset prices S , which from now onwe refer to as the unaffected asset prices, depending on order size and market conditions. These impactfunctions have been well-studied in the empirical literature, see, e.g., Easley and O’hara [1987], Kyle[1985] and the theoretical literature, see, e.g., Bertsimas and Lo [1998]. We will assume the price impactfunction from Almgren [2012], extended to multiple assets. In this framework the perceived asset prices,denoted by ˜ S , are ˜ S ( t ) = S ( t ) + Ξ( t ) v ( t ) , (2) h eq:perc˙asset i where Ξ( t ) = η ( t ) η ( t ) · · · η n ( t ) η ( t ) η ( t ) · · · η n ( t )... ... . . . ... η n ( t ) η n ( t ) · · · η nn ( t ) (3)is the matrix containing the (stochastic) instantaneous market impact coefficients η ij ( t ) ≥
0. We willassume that the diagonal is strictly positive, η ii ( t ) >
0. In our framework we assume that trading hasno impact on the market impact coefficients, i.e., η ij ( t ) moves independently of v . We want Ξ( t ) to bepositive definite for reasons explained below, which means it is a symmetric matrix. Under this modela higher impact for a given trading rate corresponds to a less liquid asset, as the trade will eat upmore of the order book (which is then instantaneously replenished). On the other hand, a lower impactcorresponds to a more liquid asset. Therefore, we will call the η liquidity parameters, which is moreconcise.In this paper we will assume that both the n volatility and n ( n + 1) / n ( n + 3) / dξ k ( t ) = − ξ k ( t ) δ k dt + β k √ δ k dB k ( t ) , j = 1 , . . . , n ( n + 3)2 . Here δ k is the market relaxation time and β k describes the dispersion of volatility ( k = 1 , . . . , n ) andliquidity ( k = n + 1 , . . . , n ( n + 3) /
2) around their average levels. The Brownian motions are correlatedas dB k ( t ) dB m ( t ) = ̺ km dt , and there is no correlation between the processes B and W . The volatilitiesdepend on these processes as σ k ( t ) = ¯ σ k e ξ k ( t ) , k = 1 , . . . , n, σ k is the average level of the k th volatility. The liquidity parameters depend similarly on theprocesses as η kℓ ( t ) = ¯ η kℓ e ξ n + k +( ℓ − n ( t ) , k = 1 , . . . , n, ℓ = 1 , . . . , k, where ¯ η kℓ is the average level of the liquidity process η kℓ . The choice of the driving stochastic processesis motivated by Almgren [2012], however, all our results are still valid for a different choice of drivingprocesses. Coordinated variation
In order to reduce the dimensionality, Almgren [2012] (which only considers n = 1) assumes that η ( t )and σ ( t ) vary perfectly inversely, i.e., σ ( t ) η ( t ) = ¯ σ ¯ η, or equivalently, ξ ( t ) = − ξ ( t ) . In terms of the dynamics of the processes, coordinated variation corresponds to setting δ = δ , β − β =0, ρ = − ξ (0) = ξ (0). This relationship is argued to be a natural consequence of a trading timemodel in which the single source of uncertainty is the arrival rate of trade events. In such a model eachtrade event brings a fixed amount of price variance and the opportunity to trade a fixed number of sharesfor a particular cost simultaneously. It is argued in the same paper that this assumption can be seriouslyviolated during events when volatility sharply increases while liquidity is withdrawn simultaneously. Wewill therefore not make this assumption. h sec:CoT i The cost of trading given a control v , denoted by C , is the difference between the amount paid to tradethe assets and its initial market value x T S (0). By using partial integration for c`adl`ag processes we find C = Z T v T ( t ) ˜ S ( t ) dt − x T S (0)= n X k =1 Z T σ k ( t ) x k ( t ) dW k ( t ) + Z T v T ( t )Ξ( t ) v ( t ) dt. We will determine the optimal control v by the mean-variance criterionmin v [ E ( C ) + λ Var( C )] , (4) h MVCriterion i where λ ≥ λ = 0 corresponds to a risk-neutral trader. Weneed to calculate the variance term in (4). Strictly speaking, it involves contributions from W k ( t ), theuncertainty in the asset price, as well as from Ξ( t ) and Σ( t ), the uncertainties in the market condition.We can circumvent the need to approximate the contributions of Ξ( t ) and Σ( t ) by making the so-calledsmall-impact approximation (see Almgren [2012]): the variance comes primarily from the price volatilityrepresented by Σ, with lesser contributions from the uncertainty in Ξ( t ) and v ( t ),Var( C ) ≈ E Z T x T ( t )Σ( t ) x ( t ) dt. This is true if the portfolio is small enough such that price changes due to impact of trading are smallcompared to volatility. Under this condition the mean-variance cost function for a control v ( t ) is E ( C ) + λ Var( C ) ≈ E Z T (cid:2) v T ( t )Ξ( t ) v ( t ) + λ x T ( t )Σ( t ) x ( t ) (cid:3) dt. (5) h eq:CCcost i v such that the above cost functionis minimized. In general, starting at some time t ≥ x ( t ) assets left to trade, we take as the valuefunction c ( t, x, Ξ , Σ) = min v ( s ) , t ≤ s ≤ T E Z Tt (cid:2) v T ( s )Ξ( s ) v ( s ) + λ x T ( s )Σ( s ) x ( s ) (cid:3) ds. (6) h eq:OptimizProblem i Since both Ξ( s ) and Σ( s ) are positive definite, the cost of trading will be positive. h sec:CC i In this section, we will only consider the static problem, i.e., the case when Ξ( t ) and Σ( t ) are constantover time. We will first show that this problem has a unique minimizer. Proposition 1.
Assume that Ξ and Σ are two positive definite matrices. Then the optimization problem min v ( t ) , ≤ t ≤ T Z T L ( t, x , v ) dt where L ( t, x , v ) = v T ( t )Ξ v ( t ) + λ x T ( t )Σ x ( t ) , has a unique minimizer in W , ([0 , T ] , R n ) Proof.
It is clear that L ( t, x , v ) is convex in v , since v T ( t )Ξ v ( t ) is a quadratic form with Ξ positivedefinite. Furthermore, since Ξ and Σ are two positive definite matrices L ( t, x , v ) ≥ α v T ( t ) v ( t ) + λα x T ( t ) x ( t ) , where α and α are the smallest eigenvalues of Ξ and Σ respectively. q.e.d. The Euler–Lagrange equations are ∂ L ∂x k − ddt ∂ L ∂x ′ k = 0 , k = 1 , . . . , n, where the total derivative equals ddt ∂ L ∂x ′ k = ∂∂t ∂ L ∂x ′ k + n X ℓ =1 x ′ ℓ ∂∂x ℓ ∂ L ∂x ′ k + n X ℓ =1 x ′′ ℓ ∂∂x ′ ℓ ∂ L ∂x ′ k , k = 1 , . . . , n. Straightforward calculations then show that the Euler–Lagrange equations reduce to the systemΞ x ′′ ( t ) = λ Σ x ( t ) . (7) h eq:EL˙eqns ih prop:CC i Proposition 2.
Under the assumption that Ξ( t ) = Ξ and Σ( t ) = Σ are constant for ≤ t ≤ T , theboundary value problem x ′′ ( t ) = λ Ξ − Σ x ( t ) , for ≤ t ≤ T, where x (0) = x and x ( T ) = , (8) h eq:sysODE i has the solution, for λ > , x ( t ) = Ω( t, T, Ξ , Σ) x , where Ω( t, T, Ξ , Σ) = sinh( C ( T − t )) sinh( CT ) − , (9) h eq:xCC i and v ( t ) = Ω ′ ( t, T, Ξ , Σ) x , where Ω ′ ( t, T, Ξ , Σ) = − cosh( C ( T − t )) sinh( CT ) − C, (10) h eq:vCC i where C is a matrix square root such that C = λ Ξ − Σ and sinh and cosh are matrix functions, i.e., sinh( A ) = P k ≥ A k +1 / (2 k + 1)! and cosh( A ) = P k ≥ A k / (2 k )! , and sinh( CT ) − is the matrix inverseof sinh( CT ) . roof. Set B = λ Ξ − Σ. This matrix has full rank by assumption and thus we can find a square root C = B / such that C = B .Now introduce y ( t ) = ( x ( t ) , x ′ ( t )) then we obtain a first order system with constant coefficients y ′ ( t ) = (cid:18) I n B (cid:19)| {z } =: A y ( t ) = A y ( t ) (11) h eq:dy i for which the fundamental matrix of solutions is given by the matrix exponentialΨ( t ) = e At = X k ≥ A k t k k ! . We have A k = B k/ B k/ ! for k = 0 , , , . . . , B ( k − / B ( k +1) / ! for k = 1 , , , . . . . Now using cosh( A ) = P k ≥ A k / (2 k )! and sinh( A ) = P k ≥ A k +1 / (2 k + 1)!, we can write, with C = B / , Ψ( t ) = cosh( Ct ) sinh( Ct ) C − sinh( Ct ) C cosh( Ct ) ! , and the general solution to (11) is thus given by y ( t ) = Ψ( t ) c . Write c = ( c , c ) then the boundary conditions of the original problem (8) give ( x (0) = cosh( C c + sinh( C C − c = x x ( T ) = cosh( CT ) c + sinh( CT ) C − c = from which it immediately follows that c = x and c = − C sinh( CT ) − cosh( CT ) x . Here sinh( CT ) − means the matrix inverse of sinh( CT ). The solution to (8) is thus given by x ( t ) = (cid:0) cosh( Ct ) − sinh( Ct )(sinh( CT )) − cosh( CT ) (cid:1) x = sinh( C ( T − t )) sinh( CT ) − x , where we used the fact that sinh( CT ) − commutes with cosh( CT ) assuming C has full rank, which followsfrom B having full rank. That the two matrices commute follows easily by using the Taylor series forsinh( A ) and cosh( A ) and the eigenvalue decomposition of A = CT = V Λ V − , i.e., cosh( A ) sinh( A ) − = V cosh(Λ) V − V sinh(Λ) − V − and the diagonal matrices sinh(Λ) − and cosh(Λ) commute. q.e.d. The above result was formulated for the interval 0 ≤ t ≤ T . The obvious change for t ≤ s ≤ T gives,with x ( t ) = x t given, x ( s ) = Ω( s − t, T − t, Ξ , Σ) x t , and v ( s ) = Ω ′ ( s − t, T − t, Ξ , Σ) x t , (12) h eq:CC i where as Ξ( s ) and Ω( s ) are constant we might think of Ξ = Ξ( t ) and Σ = Σ( t ). This will become of usein the following.These strategies do not adapt to fluctuations in market parameters, and we will denote them by CC,short for constant coefficients . For the cases n = 1 and n = 2, we have the following two corollaries.6 cor:1A i Corollary 1.
Assuming n = 1 , i.e., there is only trading in one asset, and η ( t ) = η and σ ( t ) = σ constant, the optimal trading trajectory for ≤ t ≤ T is given by x ( t ) = sinh( µ ( T − t ))sinh( µT ) x , and v ( t ) = − cosh( µ ( T − t ))sinh( µT ) µ x , with µ = λσ /η . h cor:2A i Corollary 2.
When trading in two assets, i.e., n = 2 , and Ξ( t ) = Ξ and Σ( t ) = Σ constant, the optimaltrading trajectories for ≤ t ≤ T are given by x ( t ) = 1 θ − θ θ s ( t ) − θ s ( t ) θ θ ( s ( t ) − s ( t )) s ( t ) − s ( t ) θ s ( t ) − θ s ( t ) ! x , (13) h eq:optvCC2A i and v ( t ) = 1 θ − θ θ s ′ ( t ) − θ s ′ ( t ) θ θ ( s ′ ( t ) − s ′ ( t )) s ′ ( t ) − s ′ ( t ) θ s ′ ( t ) − θ s ′ ( t ) ! x , where for λ Ξ − Σ = (cid:18) a bc d (cid:19) = 2 λη + 2 η η + η − η η (cid:18) σ σ ρ ( η + η ) − σ η σ ( η + η ) − ρσ σ η σ ( η + η ) − ρσ σ η σ σ ρ ( η + η ) − σ η (cid:19) we set D = a + 4 bc − ad + d and α = a − d, β = a + d,µ = β − √ D , µ = β + √ D ,θ = α − √ D c , θ = α + √ D c ,s ( t ) = sinh( µ ( T − t ))sinh( µ T ) , s ( t ) = sinh( µ ( T − t ))sinh( µ T ) ,s ′ ( t ) = − µ cosh( µ ( T − t ))sinh( µ T ) , s ′ ( t ) = − µ cosh( µ ( T − t ))sinh( µ T ) . Proof.
The proof follows by making an eigenvalue decomposition of B = λ Ξ − Σ = V Λ V − , then C = V Λ / V − and sinh( C ( T − t )) = V sinh(Λ / ( T − t )) V − and sinh( CT ) − = V sinh(Λ / T ) − V − . For2 × q.e.d. Finally, we have the following corollary. h cor:CC i Corollary 3.
Assume that Ξ( s ) = Ξ and Σ( s ) = Σ are constant over t ≤ s ≤ T . Denote with x ( s ) the optimal trading strategy for t ≤ s ≤ T found using (7) , assuming the asset level at time s = t to be x t . As a consequence of Proposition , the cost function (5) over [ t, T ] is a quadratic polynomial in theelements of x t .Proof. Due to Proposition 2, the cost function for this strategy, with Ω( s ) = Ω( s − t, T − t, Ξ , Σ) andΩ ′ ( s ) = Ω ′ ( s − t, T − t, Ξ , Σ), is given by Z Tt (cid:2) v T ( s )Ξ v ( s ) + λ x T ( s )Σ x ( s ) (cid:3) ds = x Tt Z Tt (cid:2) Ω ′ ( s ) T ΞΩ ′ ( s ) + λ Ω( s ) T ΣΩ( s ) (cid:3) ds ! x t = x Tt Q x t and is a sum of quadratic expressions in x t with coefficients that are integrals independent of x t . There-fore the cost function when assuming constant Ξ and Σ is a quadratic polynomial in the initial assetposition x ( t ). q.e.d. The dynamic problem h sec:DynProb i This section deals with the more general situation where both Ξ( t ) and Σ( t ) move stochastically overtime. In this case, the Bellman principle can be used on the value function (6), c ( t, x, Ξ , Σ) = min v ( t ) (cid:2) v T ( t )Ξ( t ) v ( t ) dt + λ x T ( t )Σ( t ) x ( t ) dt + E c ( t + dt, x + dx, Ξ + d Ξ , Σ + d Σ) (cid:3) . (14) h eq:optprob i Unfortunately finding an analytic solution to the problem is impossible in a general setting, and numericaltechniques have to be used. The case of one asset under coordinated variation, studied in Almgren [2012],reduces to a PDE with one spatial dimension. Solving this problem can be done efficiently using stan-dard techniques, and was done using a finite difference scheme in the aforementioned paper. Increasingthe number of assets and relinquishing the coordinated variation condition quickly results in a multidi-mensional problem that is impractical for finite difference techniques (see also Longstaff and Schwartz[2001]). h sec:RHA i The so-called rolling horizon strategy (RHS) proposed in Almgren [2012] offers a dynamic, but suboptimaltrading strategy that does not need any numerical algorithm to compute. The idea is to plug in theinstantaneous values of of Ξ( t ) and Σ( t ) in the static solution, i.e., the solution to (7). From (12) we canwrite the instantaneous trading rate under RHS as v ( t ) = Ω ′ (0 , T − t, Ξ( t ) , Σ( t )) x ( t ). The algorithmtherefore assumes that the current market parameters will remain constant over the remainder of theprogram. When they change, the trading speed is altered using the new values. The author argues thatthis strategy is strictly optimal only in the infinite-horizon case, and only when the market parameterscovary in the appropriate way (but without specifying how). It is furthermore claimed to provide areasonable approximation, that is easy to implement.We start by proving the following proposition which states under what condition the RHS reduces tothe static solution. Proposition 3.
Consider the execution problem for n assets with initial position x . If λ Ξ − Σ → n ,with n the n × n zero matrix, then the RHS converges to the CC solution. h prop:RHS2 i Proof.
When λ Ξ − Σ → n , the system of ODEs (8) reduces to x ′′ ( t ) = 0 , which leads to the solution x ( t ) = x (cid:18) − tT (cid:19) , i.e., the trading happens linearly independent of even the current market parameters. Consequently, theRHS coincides with the static solution. q.e.d. Conversely, the following proposition (which to our knowledge has not been studied in the literature)shows when the RHS becomes optimal, at least for the case n = 1. Proposition 4.
Assume n = 1 . Then for λ ¯ σ / ¯ η → ∞ , the rolling horizon approximation converges tothe optimal solution if and only if ̺β β r δ δ δ − ξ + 12 δ δ β − δ δ ξ + 18 β δ = 0 . (15) h eq:RHAcond i Note that this equation is satisfied under coordinated variation. h prop:RHS i Proof.
The optimization problem (6) for the case n = 1 can be written as c ( t, x, ξ , ξ ) = min v ( s ) , t ≤ s ≤ T E Z Tt (cid:2) η ( s ) v ( s ) + λσ ( s ) x ( s ) (cid:3) ds. c t − ξ δ c ξ + β δ c ξ ξ − ξ δ c ξ + β δ c ξ ξ + ̺ β β √ δ δ c ξ ξ + λ ¯ σ e ξ x + min v h vc x + ¯ ηe ξ v i = 0 . The minimum is v = − c x e − ξ η , which means the PDE for c is c t − ξ δ c ξ + β δ c ξ ξ − ξ δ c ξ + β δ c ξ ξ + ̺ β β √ δ δ c ξ ξ + λ ¯ σ e ξ x − c x e − ξ η = 0 . The value function c is strictly proportional to x , which means we can nondimensionalize using δ asthe time scale and τ = ( T − t ) /δ , c ( t, x, ξ , ξ ) = ¯ ηx δ u ( τ, ξ , ξ ) , where u is a nondimensional function of nondimensional variables. The PDE becomes u τ + ξ u ξ + δ δ ξ u ξ − ¯ µ δ + e − ξ u − β u ξ ξ − ̺ r δ δ β β u ξ ξ − β u ξ ξ = 0 , where ¯ µ = λ ¯ σ / ¯ η . We assume that this PDE has a unique solution. The trade velocity in terms of thetransformed value function is v = xδ e − ξ u ( τ, ξ , ξ ) . Using Corollary 1 the continuous time rolling horizon strategy is given as − x ¯ µe ξ − ξ coth (cid:16) ¯ µe ξ − ξ ( T − t ) (cid:17) , with ¯ µ = λ ¯ σ / ¯ η . In terms of u this gives u = − δ ¯ µe ξ − ξ coth (cid:16) ¯ µe ξ − ξ ( T − t ) (cid:17) . Filling this in in the PDE for u gives the equation (cid:18) δ δ β − β δ + ̺β β p δ δ δ (cid:19) τ e ξ − ξ coth (cid:16) ¯ µe ξ − ξ δ τ (cid:17) (cid:16) − coth (cid:16) ¯ µe ξ − ξ δ τ (cid:17)(cid:17) ¯ µ + − δ δ ξ τ − δ + 12 ξ δ τ − β δ τ − δ + 3 δ δ β − ̺β β s δ δ δ τ ! e ξ (cid:16) − coth (cid:16) ¯ µe ξ − ξ δ τ (cid:17)(cid:17) ¯ µ + ̺β β s δ δ δ − ξ + 12 δ δ β − δ δ ξ + 18 β δ ! e ξ + ξ coth (cid:16) ¯ µe ξ − ξ δ τ (cid:17) ¯ µ. It is clear that in order for this equation to hold when λ ¯ σ / ¯ η → ∞ , we need12 ̺β β r δ δ δ − ξ + 12 δ δ β − δ δ ξ + 18 β δ = 0 , which is exactly (15). Under coordinated variation, δ = δ , β − β = 0, ̺ = − ξ + 2 ξ = 0,which satisfies the above equation. q.e.d. Numerical experiments, some of which will be shown in Section 5.3, implicate that λ ¯ σ / ¯ η does noteven need to be that large for the rolling horizon solution to become (nearly) optimal. It is not clearwhat happens in the multi-asset case. The results in Section 5 seem to indicate that the RHS can still9ecome optimal for certain parameter choices, but finding the exact condition lies outside the scope ofthis paper.In the remainder of this section we will discuss the RHS in a discretized time framework, which willprove useful in the next section. Suppose we discretize time as ( t k ) k =0 ,...,M , using step length ∆ t . Fromnow on, the notation x k = x ( t k ) and v k = v ( t k ) will be used. Assume we are at time t k with asset levels x k and Ξ( t k ) and Σ( t k ) known (i.e., observed). We need to decide on the trading speed over the interval[ t k , t k +1 ] which will be taken constant and is denoted by v k . To decide on the optimal value of v k wewill assume that Ξ( t k ) and Σ( t k ) remain constant over the remainder of the program [ t k , T ] and use theCC solution. In fact, by using (9) we know the optimal value of x k +1 directly, x k +1 = Ω(0 , T − t k , Σ( t k ) , Ξ( t k )) x k = Ω k x k where the shorthand notation Ω k = Ω(0 , T − t k , Σ( t k ) , Ξ( t k ))was introduced. The trading rate can be deduced as v k = x k +1 − x k ∆ t = (Ω k − I n ) x k ∆ t . The trader will repeat this process at the next time step.Each time step costs calculating Ω k and a matrix-vector product. If we assume the cost of calculatingΩ k to be O ( n ), e.g., by making use of an eigen decomposition, then this dominates the cost per step.The cost per step is thus O ( n ).If we assume we know all values of Ξ( t k + ℓ ) and Σ( t k + ℓ ) then we could propagate this as an iterativescheme and write, for ℓ > x k + ℓ = Ω k + ℓ − · · · Ω k x k = V ( k, ℓ, Ξ( · ) , Σ( · )) Ω k x k = V ( k, ℓ, Ξ( · ) , Σ( · )) x k +1 , (16) h eq:xkl i where the propagation from x k +1 to x k + ℓ , ℓ ≥
1, is given by V ( k, ℓ, Ξ( · ) , Σ( · )) = ℓ − Y m =1 Ω k + ℓ − m . The reason to define a propagation matrix from x k +1 to x k + ℓ instead of from x k to x k + ℓ has to do withthe RHMC scheme which will be described next. We note that V ( k, ℓ + 1 , Ξ( · ) , Σ( · )) = Ω k + ℓ V ( k, ℓ, Ξ( · ) , Σ( · )) . (17) h eq:Vrecur i An overview of the discrete RHS method is given in Algorithm 1.The following proposition ensures the RHS scheme is numerically stable.
Proposition 5.
The RHS scheme (16) is numerically stable, i.e., all the eigenvalues of Ω k = Ω(0 , T − t k , Σ( t k ) , Ξ( t k )) are positive and bounded by one for all t k ∈ [0 , T ] . h prop:stability i Proof.
From Proposition 2 we know thatΩ k = Ω(0 , T − t k , Σ( t k ) , Ξ( t k ))= sinh( C ( T − t )) sinh( CT ) − where C is a square root of B = λ Ξ( t k ) − Σ( t k ). The dependence of C and B on t k is suppressedfor ease of notation. The proof follows by making an eigenvalue decomposition of B = V Λ V − , then10 = V Λ / V − and sinh( C ( T − t )) = V sinh(Λ / ( T − t )) V − and sinh( CT ) − = V sinh(Λ / T ) − V − .Therefore Ω k = V sinh(Λ / ( T − t )) sinh(Λ / T ) − V − , which implies that the eigenvalues of Ω k are positive and bounded by one, since sinh is a positive andstrictly increasing function on (0 , ∞ ). Therefore the eigenvalues of the matrix product V in (17) arealso positive and bounded by one. Because the time step k was arbitrary, the stability holds for all timesteps, leading to the conclusion that the RHS algorithm is stable over [0 , T ]. q.e.d.Algorithm 1 The discrete RHS for optimal asset execution. h alg:RHS i Fix a set of M + 1 trading times, equally spaced with length ∆ t = T /M .At time t k observe Ξ( t k ) and Σ( t k ).Calculate x k +1 = Ω k x k .Trade x k +1 − x k assets. h sec:MC i We are now ready to introduce our new method: the rolling horizon Monte Carlo (RHMC) method. Wediscretize time again as ( t k ) k =0 ,...,M , using step length ∆ t , with x k = x ( t k ) and v k = v ( t k ). Using theBellman principle the optimization problem (14) at time step t k becomesmin v k (cid:26) v Tk Ξ( t k ) v k ∆ t + λ x Tk Σ( t k ) x k ∆ t + E [ c ( x k +1 , Ξ( t k +1 ) , Σ( t k +1 )) | x k , v k , Ξ( t k ) , Σ( t k ) ] (cid:27) , where c is the total cost incurred from time t k +1 onwards, assuming the control is used. We firstexplain the standard technique to solve this problem in a backwards manner. Since the expectation isconditional on the current values of x k , v k , Ξ( t k ) and Σ( t k ) one has to account for all possible valuesof these conditional parameters. When using Monte Carlo, one samples Ξ( t ) and Σ( t ) at time steps t , . . . , t M . The continuation value , i.e., the expected value of c in the optimization problem, is thenapproximated using a multivariate regression on powers or exponentials of the underlying variables.However, one has to be careful since the variable x k is endogenous : it is determined completely by thecontrol v k . As in Boogert and De Jong [2008] we notice that the continuation value depends only onthe asset level that is reached, not on the previous level and chosen trading rate. Therefore, one coulddiscretize the endogenous variable x and do separate regressions for each level, depending only on theprocesses ξ . This scheme still needs an exponentially increasing number of regressions as the number ofassets increases. We therefore propose to use a rolling horizon Monte Carlo scheme that we will showdoes not require the discretization of the endogenous variable, nor the use of regressions.While the previous methods operate backwards , the rolling horizon Monte Carlo algorithm we proposeis a forward scheme that approximates the continuation value using a sub-optimal control. As a majoradvantage we can now assume that the conditional values x k , v k , Ξ( t k ) and Σ( t k ) are known and we areonly left with the evaluation of a single expectation. RHMC-I: Rolling horizon Monte Carlo with RHS until the end
Assume we are at time t k with asset levels x k and Ξ( t k ) and Σ( t k ) known (i.e., observed). We need todecide on the trading speed over the interval [ t k , t k +1 ] which will be taken constant and is denoted by11 k . We will try to find the optimal value of v k , denoted by v ∗ k , that minimizes the cost c ( t k , x k , Ξ( t k ) , Σ( t k ))= min v ( s ) t k ≤ s In each time step t k of the discretization, given the sampled matrix A k +1 , there existsa unique optimal position x k +1 , which is given as the solution to (cid:0) Ξ( t k ) + A k +1 (∆ t ) (cid:1) x k +1 = Ξ( t k ) x k . (21) h eq:optControl ih prop:optimalv i Proof. We reduce the notational overload in the minimization expression by setting x = x k +1 , Ξ = Ξ( t k ), A = A k +1 (Ξ( t k ) , Σ( t k )). We need to minimize the function f ( x ) = x T (cid:18) Ξ∆ t + A ∆ t (cid:19) x − ( x T Ξ x k + x Tk Ξ x ) 1∆ t . Setting the gradient with respect to x to zero we obtain ∇ f ( x ) = 2 (cid:18) Ξ∆ t + A ∆ t (cid:19) x − t x k = , and thus x is the solution to (cid:0) Ξ + A (∆ t ) (cid:1) x = Ξ x k , provided there exists a unique minimum. This can be checked by the positive definiteness of the Hessianmatrix which is given by H ( f ) = 2 (cid:18) Ξ ∆ t + A ∆ t (cid:19) . Remember that Ξ is positive definite by assumption. We next show that A = A k +1 is positive definite.We have that A k +1 is the sum of matrices Q k,ℓ (Ξ( · ) , Σ( · )) = V T ( k, ℓ, Ξ( · ) , Σ( · )) " Ω Tk + ℓ − I n ∆ t Ξ( t k + ℓ ) Ω k + ℓ − I n ∆ t + λ Σ( t k + ℓ ) V ( k, ℓ, Ξ( · ) , Σ( · )) . For ℓ = 1 we have V ( k, ℓ, Ξ( · ) , Σ( · )) = I n and thus Q k, (Ξ( · ) , Σ( · )) = Ω Tk +1 − I n ∆ t Ξ( t k +1 ) Ω k +1 − I n ∆ t + λ Σ( t k +1 ) . 13y Proposition 5 we know that the matrices Ω k +1 have eigenvalues in the interval (0 , k +1 − I n has eigenvalues in the interval ( − , Tk +1 − I n ) Ξ( t k +1 ) (Ω k +1 − I n ) positive definite, since Ξ( t k +1 ) is positive definite by construction. For ℓ > k + ℓ − m to the right and to the left. Therefore A k +1 is positive definite. q.e.d.Proposition 7. The RHMC-I algorithm is stable, i.e., the eigenvalues of (cid:0) Ξ( t k ) + A k +1 (∆ t ) (cid:1) − Ξ( t k ) (22) h eq:stableRHMC i are positive and bounded by one for all k . h prop:stableRHMC i Proof. Denote by µ , . . . , µ n the eigenvalues of (22) in increasing order. Because the product of twopositive definite matrices has positive eigenvalues (see Horn and Johnson [1985]), we have µ > µ n (cid:16)(cid:0) Ξ( t k ) + A k +1 (∆ t ) (cid:1) − Ξ( t k ) (cid:17) ≤ µ n (cid:16)(cid:0) Ξ( t k ) + A k +1 (∆ t ) (cid:1) − (cid:17) µ n (Ξ( t k ))= (cid:0) µ (cid:0) Ξ( t k ) + A k +1 (∆ t ) (cid:1)(cid:1) − µ n (Ξ( t k )) . From the Weyl inequalities Bhatia [2001] we have that (cid:0) µ (cid:0) Ξ( t k ) + A k +1 (∆ t ) (cid:1)(cid:1) − ≤ (cid:0) µ n (Ξ( t k )) + (∆ t ) µ ( A k +1 ) (cid:1) − , leading to the inequality µ n (cid:16)(cid:0) Ξ( t k ) + A k +1 (∆ t ) (cid:1) − Ξ( t k ) (cid:17) ≤ µ n (Ξ( t k )) µ n (Ξ( t k )) + (∆ t ) µ ( A k +1 ) ≤ A k +1 , which was proven in Proposition 6. q.e.d. RHMC-II: Rolling horizon Monte Carlo with CC till the end Under the RHMC-I method the trading rate in step t k was decided by using the RHS scheme andcalculating the matrix A k +1 (Ξ( · ) , Σ( · )) , and then using Proposition 6. It is of course not mandatoryto use RHS as the suboptimal control. We could just as well use the CC to determine A k +1 (Ξ( · ) , Σ( · ))which will reduce the computational complexity. The difference with RHMC-I therefore is that we nowassume the matrices Ξ( s ) and Σ( s ) to be fixed from time t k +1 on to calculate the suboptimal control.We can therefore use (12) to write (19), for 1 ≤ ℓ < M − k , as v Tk + ℓ Ξ( t k + ℓ ) v k + ℓ + λ x Tk + ℓ Σ( t k + ℓ ) x k + ℓ = x Tk +1 Q k,ℓ (Σ( · ) , Ξ( · )) x k +1 , where Q k,ℓ is now defined as Q k,ℓ (Σ( · ) , Ξ( · )) = Ω ′ Tk,ℓ Ξ( t k + ℓ ) Ω ′ k,ℓ + λ Ω Tk,ℓ Σ( t k + ℓ ) Ω k,ℓ (23) h eq:Q˙rhmc2 i and Ω k,ℓ = Ω( t k + ℓ − t k +1 , T − t k +1 , Ξ( t k +1 ) , Σ( t k +1 )). Both Proposition 6 and Proposition 7 are stillvalid for the RHMC-II method.Assuming the left rectangle rule we arrive at a cost O ( N ( M − k ) n ) which is the same complexity asRHMC-I but it will in practice be (much) faster since the matrix Ω can now be calculated as a function oftime instead of being updated on each of the future trading dates and there is no more need to computethe function V which is essentially a cumulative matrix product.An overview of both of the algorithms is given in Algorithm 2. Note that in the single asset case undercoordinated variation and for λ ¯ σ / ¯ η → ∞ our RHMC-I scheme converges to the optimal solution sinceProposition 4 shows that the RHS converges to the optimal solution, and the continuation values arecalculated using this optimal control. Other than this case we expect RHMC-I to improve significantlyon the rolling horizon solution for any choice of market parameters. These outlooks can be justified sincethe RHS only considers the current market conditions and gives the trading rate under the assumption14hat these conditions remain constant for the remaining trading period. Our algorithm, on the otherhand, uses a Monte Carlo procedure to gain information about future trading, and chooses a tradingrate accordingly.It is at this point not clear whether the RHMC-II method will provide an improvement over theRHS. At least in the case n = 1 under coordinated variation and λ ¯ σ / ¯ η → ∞ we expect the RHMC-IImethod to be worse than the RHS, since the RHS will tend to be optimal whereas the CC will not.The numerical results in Section 5 show that the RHMC-II algorithm indeed underperforms in this case,but as soon as coordinated variation is not assumed, it performs similarly to the RHMC-I algorithm.Interestingly, in the case n = 2 it appears that there is almost no difference between the RHMC-I andRHMC-II algorithms for all market parameters considered. Algorithm 2 The Monte Carlo rolling horizon method for optimal asset execution. h alg:MCRHalg i Fix a set of M + 1 trading times, equally spaced with length ∆ t = T /M .At time t k observe Ξ( t k ) and Σ( t k ).Sample N paths for Ξ( · ) and Σ( · ) given Σ( t k ) and Ξ( t k ) using (quasi-)Monte Carlo.Calculate the matrix A (20) by using (18) for the RHMC-I method or (23) for the RHMC-II method.Calculate x k +1 using (21).Trade x k +1 − x k assets. h sec:OPT i It is possible to compute the optimal control for the dynamic problem given a discretization and assumingthe paths of the market parameters are known. This is of course of no use for a trader facing an executionproblem, but it allows us to check how the cost of trading is situated for all methods against the costwhen trading with the optimal discrete control.Assume for now that there is only trading in one asset, i.e., n = 1. We use again the discretizationof time ( t k ) k =0 ,...,M with step length ∆ t = T /M , and approximate the integrals using a left-point rule.Given the paths of the market parameters, Ξ( t k ) and Σ( t k ), k = 0 , . . . , M , the optimal control problemcorresponds to solving the minimization problemminimize v k ,k =0 ,...,M ∆ t M − X k =0 v k Ξ( t k ) + λx k Σ( t k )subject to x = X, By noting that x ( t ) = − Z Tt v ( s ) ds, we find for the discretized points x ( t k ) x k = − M − X m = k v m ∆ t. Therefore, the minimization problem can be written asminimize v k ,k =0 ,...,M ∆ t M − X k =0 v k Ξ( t k ) + λ ∆ t Σ( t k ) M − X m = k v m ! subject to − M − X m =0 v m ∆ t = X, v k ,k =0 ,...,M ∆ t M − X k =0 v k Ξ( t k ) + λ ∆ t M − X k =0 M − X ℓ =0 v k v ℓ min( k,ℓ ) X m =0 Σ( t m )subject to − M − X m =0 v m ∆ t = X, By introducing the matrix ˜Σ where ˜Σ kℓ = P min( k,ℓ ) m =0 Σ( t m ), the diagonal matrix ˜Ξ where ˜Ξ kk = Ξ( t k )and the vector v t which is the time-discretized control v , i.e., v tm = v m , we can rewrite the above problemto minimize v t ( v t ) T (cid:16) ∆ t ˜Ξ + λ ∆ t ˜Σ (cid:17) v t subject to − ∆ t TM v t = X, Since the matrix ∆ t ˜Ξ + λ ∆ t ˜Σ is positive definite, a unique minimizer exists.The more general case n > v tk the time-discretizedvectors of asset k , k = 1 , . . . , n The minimization problem becomesminimize v t ,...,n ( v t ,...,n ) T ∆ t ˜Ξ (11) · · · ˜Ξ (1 n ) ... . . . ...˜Ξ ( n · · · ˜Ξ ( nn ) + λ ∆ t ˜Σ (11) · · · ˜Σ (1 n ) ... . . . ...˜Σ ( n · · · ˜Σ ( nn ) v t ,...,n subject to − ∆ t w T v t ,...,n = X , ... − ∆ t w Tn v t ,...,n = X n , where ˜Ξ ( ij ) is the diagonal matrix with ˜Ξ ( ij ) kk = Ξ ij ( t k ), ˜Σ ( ij ) kℓ = P min( k,ℓ ) m =0 Σ ij ( t m ), v t ,...,n = v t ... v tn and w k = M ( k − M M ( n − k − . h sec:NumResults i In this section we numerically illustrate RHMC-I and RHMC-II, and compare it to the CC, RHS anddiscrete optimal solutions. The aim of our method is to outperform the RHS solution. The results inthis section show that we are successful at this goal.To choose the number of simulations used for calculating the expected continuation value, we havelooked at the cost of trading for different parameter choices in function of the number of samples used.A typical example is given in Figure 1a for n = 1 and Figure 1b for n = 2, using the RHMC-I method.It is clear that using quasi-Monte Carlo (the red curves) is advantageous over using plain Monte Carlo(the blue curves). For our methods we use a Sobol’ sequence with parameters from Joe and Kuo [2008].We will choose N = 500 for our numerical experiments, but we would also like to point out that taking N = 200 using quasi-Monte Carlo seems sufficient for practical applications. For the one asset case, assuming coordinated variation holds, we ran experiments using fixed parameters T = 10, ∆ t = 1 / β = δ = 1, and different values of ¯ σ , ¯ η and λ . The dimension of the problem, atsome arbitrary time step k ∈ { , , . . . , } is k , since there is only one stochastic process driving themarket parameters.Table 1 gives an overview of the cost of trading for different parameter values using 200 simulationruns, i.e., 200 sample paths of ξ on which we used the different execution algorithms. Since this is also16 fig:1A˙conv i , . . . . . . N C o s t o f t r a d i n g (a) One Asset h fig:2A˙conv i , . . . . . N (b) Two Assets Figure 1: The cost of trading given specific paths for the liquidity and volatility parameters, in function of N , the number of (quasi-)Monte Carlo paths used to determine the mean continuation value. Each bluecurve corresponds to a different seed for the Monte Carlo numbers used, and each red curve correspondsto a different digital shift used to obtain the quasi-Monte Carlo sample. h fig:1A2A˙conv i the setting considered in Almgren [2012], we have implemented the finite difference scheme as outlinedin said paper. We will call this solution the optimal continuous solution. It should be noted that thediscrete optimal solution and the continuous optimal solution do not need to coincide.The effect of Proposition 3 is clearly visible in the table: the RHS and CC have similar costs oftrading when λ ¯ σ / ¯ η is smaller than 10 − . It is also in these cases that both RHMC-I and RHMC-IIoutperform the RHS method significantly, reducing the extra cost over the continuous optimal solutionto about a fourth of that of RHS. This reduction in cost is constant over all parameter choices for which λ ¯ σ / ¯ η is smaller than 10 − . Both our methods have similar, if not identical, reductions in cost.On the other hand, for λ ¯ σ / ¯ η larger than 10 − the effect of Proposition 4 becomes visible in the results,with the extra cost of the RHS compared to the continuous optimal solution dropping significantly. OurRHMC-I method keeps outperforming the RHS method, especially for the cases where λ ¯ σ / ¯ η is of theorder 10 − . The RHMC-II method starts to lag behind the RHMC-I method for these cases, with smallerreductions in cost. For the cases where λ ¯ σ / ¯ η is of the order 10 − the advantage of the RHMC-I methodis greatly reduced, because the extra cost of the RHS method is almost zero. The RHMC-II methodeven fails to outperform the RHS method in these cases. This should not come as a surprise, since theRHMC-I method uses the RHS as suboptimal control, so if the RHS converges to the optimal control,the RHMC-I method also converges to the optimal control by construction. Because RHMC-II uses theCC as suboptimal control, it does not show this convergence.It is noteworthy to mention the increase in cost should a trader use the CC solution when marketparameters are not constant. We ran numerical experiments using fixed parameters T = 10 and β = β = δ = δ = 1, for differentvalues of ¯ σ , ¯ η , λ and ̺ , and a time discretization ∆ t = 1 / k ∈ { , , . . . , } is 2 k , since there are two stochastic processes driving the marketparameters.Table 2 gives an overview of the cost of trading for different parameter values using 200 simulationruns, i.e., 200 sample paths of Ξ and Σ on which we used the different execution algorithms. Thepercentages indicate the increase in cost compared to the discrete optimal solution. The effects ofProposition 3 are again clearly visible in the table: for values of λ ¯ σ / ¯ η of the order 10 − or smaller theRHS and CC algorithms almost coincide. Our RHMC-I and RHMC-II method significantly outperform17he RHS method, cutting the extra cost over the discrete optimal solution to less than half. Both methodshave almost identical performance.Interestingly, even for λ ¯ σ / ¯ η of order larger than 10 − both RHMC-I and RHMC-II methods signif-icantly outperform the RHS method, as opposed to the coordinated variation case. The reduction forthe RHMC-I method is still around half, whereas the RHMC-II method is now showing reductions thatare less than half of the extra cost. h sec:2A i We ran numerical experiments using fixed parameters T = 10, ¯ σ = 1 / σ = 3 / η = 1 / η = ¯ η = 1 / β k = δ k = 1 for k = 1 , . . . , 5. The correlation matrix of the Brownian motionsdriving the market parameters is fixed as ̺ = 110 10 8 1 − − 68 10 1 − − 61 1 10 − − − − − − − − . Experiments were run with different values of x , ¯ η , ρ and λ . Time is discretized with step length∆ t = 1 / k ∈ { , , . . . , } is 5 k , sincethere are five stochastic processes driving the market parameters.Table 3 shows the costs of all methods and the improvements over the RHS of our method. It alsoshows the largest absolute value in the matrix λ Ξ − Σ, an indication of how close this matrix is to thezero matrix. A first glance at the results shows the interesting observation that there seems to be almostno difference between the RHMC-I and RHMC-II algorithm for all cases considered. The results aresimilar to the one asset case when λ Ξ − Σ is fairly close to zero, i.e., the extra cost of trading over usingthe optimal control for the RHMC-I and RHMC-II methods are about one third of those of RHS for λ Ξ − Σ smaller than the order 10 − . For λ Ξ − Σ of order 10 − the RHMC-I and RHMC-II methods stillperform similarly, still outperforming the RHS significantly in all the cases considered. This is especiallytrue in the case where opposite initial positions have to be traded: under these circumstances, costsare cut by a fifth to a sixth compared to RHS for both methods. As mentioned before, there seems tobe evidence though that there is some extension to Proposition 4 to multiple assets, but there is alsoevidence that it is not sufficient to let λ Ξ − Σ → n . Also note the huge increase in cost in some caseswhen using CC, which can be up to three times the cost compared to the optimal solution.The results for two assets are again very satisfactory (even more than the one asset case). h sec:end i In this paper we have studied the static solution of the multidimensional extension to the model inAlmgren [2012]. We also presented a new result on the rolling horizon strategy (RHS), an approximatescheme constructed in the same paper to trade in a dynamic market. Our aim was to develop a methodthat performs better than the RHS, but is still easier to compute than the optimal solution. This leadto the rolling horizon Monte Carlo methods (RHMC), that determine the optimal trading rate in eachtime step by balancing the cost of trading over the current time interval against the projected futuretrading costs. The future trading costs are determined using (quasi-)Monte Carlo and a suboptimaltrading strategy, in our case the RHS (used in RHMC-I) and the constant coefficient solution (used inRHMC-II). It was argued to perform better than the RHS since the latter only considers the currentmarket information, whereas our method uses current information as well as a prediction on futurecosts. The advantages of this scheme are that it is easy to understand and there is no need to solvea high-dimensional PDE (when using for instance a finite-difference solver) or regressions (when usingleast-squares Monte Carlo). It turns out that using our method we can reduce the extra cost of tradingover that of the RHS to as much as one sixth. Furthermore, our methods seem to be more consistent incutting trading cost compared to the RHS when there is more than one asset to be traded.18t would be very interesting to see if these results carry over to other models as well, especially whenprice impact is not instantaneous but temporary, see for instance Alfonsi et al. [2012], Gatheral et al.[2012]. Another interesting expansion would be to include directional bets in the strategy, see Almgren and Lorenz[2006] and Engle and Ferstenberg [2007]. Checking what happens if the small-impact approximation isviolated, i.e., if the variance in the cost of trading also depends on the liquidity and chosen trading ratescould also lead to interesting results. Appendix: Tables ontinuous Discrete(¯ σ, ¯ η, λ ) Optimal CC RHS RHMC-I RHMC-II Optimal λσ /η (0 . . − ) 1 . 97 2 . 58 (31%) 2 . 58 (31%) 2 . 11 (7 . . 11 (7 . . 82 5 × − (0 . . . . 00 2 . 62 (30%) 2 . 60 (29%) 2 . 14 (6 . . 14 (6 . . 86 5 × − (0 . . . 1) 4 . 60 5 . 43 (17%) 4 . 73 (2 . . 62 (0 . . 69 (2 . . 56 5 × − (0 . . − ) 2 . 96 3 . 87 (31%) 3 . 87 (31%) 3 . 17 (7 . . 17 (7 . . 73 3 . × − (0 . . . . 99 3 . 91 (30%) 3 . 89 (30%) 3 . 20 (6 . . 20 (7 . . 77 3 . × − (0 . . . 1) 5 . 81 6 . 92 (19%) 6 . 08 (4 . . 85 (0 . . 94 (2 . . 71 3 . × − (0 . . − ) 1 . 97 2 . 58 (31%) 2 . 58 (30%) 2 . 11 (7 . . 11 (7 . . 82 2 × − (0 . . . . 11 2 . 72 (29%) 2 . 65 (25%) 2 . 23 (5 . . 24 (6 . . 97 2 × − (0 . . . 1) 8 . 99 10 . . 05 (0 . . 04 (0 . . 14 (1 . . 02 2 × − (0 . . − ) 2 . 96 3 . 88 (31%) 3 . 87 (30%) 3 . 17 (7 . . 17 (7 . . 74 1 . × − (0 . . . . 10 4 . 02 (29%) 3 . 94 (27%) 3 . 29 (6 . . 30 (6 . . 88 1 . × − (0 . . . 1) 11 . . . . . . . . . . × − Table 1: The cost of trading incurred for different values of ̺ , ¯ σ , ¯ η and λ , for the one asset case ( n = 1) assuming coordinated variation. The percentagevalues denote the increase compared to the continuous time optimal solution. Smaller percentages are better. These results are based on 200 realizationsof ξ . h table:ex1ACV i iscrete(¯ σ, ¯ η, λ ) Optimal CC RHS RHMC-I RHMC-II λσ /η̺ = − . . − ) 1 . 82 2 . 58 (42%) 2 . 58 (42%) 2 . 11 (18%) 2 . 11 (18%) 10 − (0 . . . . 86 2 . 62 (41%) 2 . 61 (40%) 2 . 15 (17%) 2 . 15 (17%) 10 − (0 . . . 1) 4 . 63 5 . 57 (20%) 5 . 01 (8 . . 77 (3 . . 82 (4 . − (0 . . − ) 2 . 73 3 . 87 (42%) 3 . 87 (42%) 3 . 17 (18%) 3 . 17 (18%) 1 . × − (0 . . . . 77 3 . 91 (42%) 3 . 90 (41%) 3 . 20 (18%) 3 . 20 (18%) 1 . × − (0 . . . 1) 5 . 81 7 . 09 (22%) 6 . 44 (10%) 6 . 03 (3 . . 10 (5 . . × − (0 . . − ) 1 . 82 2 . 58 (42%) 2 . 58 (42%) 2 . 11 (18%) 2 . 11 (18%) 4 × − (0 . . . . 98 2 . 73 (38%) 2 . 68 (35%) 2 . 25 (15%) 2 . 25 (15%) 4 × − (0 . . . 1) 8 . 98 10 . . 41 (4 . . 21 (2 . . 31 (4 . × − (0 . . − ) 2 . 74 3 . 88 (42%) 3 . 87 (42%) 3 . 17 (18%) 3 . 17 (18%) 5 × − (0 . . . . 89 4 . 03 (40%) 3 . 98 (38%) 3 . 30 (16%) 3 . 31 (16%) 5 × − (0 . . . 1) 11 . . . . . . . . × − ̺ = 60%(0 . . − ) 1 . 82 2 . 58 (42%) 2 . 58 (42%) 2 . 11 (18%) 2 . 11 (18%) 10 − (0 . . . . 86 2 . 62 (41%) 2 . 62 (41%) 2 . 15 (18%) 2 . 15 (18%) 10 − (0 . . . 1) 4 . 77 5 . 64 (18%) 5 . 45 (14%) 5 . 02 (5 . . 04 (6 . − (0 . . − ) 2 . 73 3 . 87 (42%) 3 . 87 (42%) 3 . 17 (18%) 3 . 17 (18%) 1 . × − (0 . . . . 78 3 . 91 (42%) 3 . 91 (42%) 3 . 21 (18%) 3 . 21 (18%) 1 . × − (0 . . . 1) 5 . 98 7 . 18 (20%) 6 . 95 (16%) 6 . 31 (6 . . 34 (6 . . × − (0 . . − ) 1 . 82 2 . 58 (42%) 2 . 58 (42%) 2 . 11 (18%) 2 . 11 (18%) 4 × − (0 . . . . 99 2 . 74 (38%) 2 . 72 (37%) 2 . 26 (16%) 2 . 26 (16%) 4 × − (0 . . . 1) 9 . 31 10 . . . . 79 (5 . . 82 (6 . × − (0 . . − ) 2 . 74 3 . 88 (42%) 3 . 88 (42%) 3 . 17 (18%) 3 . 17 (18%) 5 × − (0 . . . . 90 4 . 03 (40%) 4 . 02 (39%) 3 . 32 (16%) 3 . 32 (17%) 5 × − (0 . . . 1) 11 . . . . . . . × − Table 2: The cost of trading incurred for different values of ̺ , ¯ σ , ¯ η and λ , for the one asset case ( n = 1). The percentage values denote the increasecompared to the optimal solution. Smaller percentages are better. These results are based on 200 realizations of ξ and ξ . h table:ex1A i iscrete(¯ σ, ¯ η, λ ) Optimal CC RHS RHMC-I RHMC-II max { λ Ξ − Σ } x = (100 , T ( − . . − ) 4 . 23 5 . 46 (29%) 5 . 46 (29%) 4 . 63 (9 . . 63 (9 . . × − ( − . . . . 28 5 . 51 (28%) 5 . 47 (28%) 4 . 69 (9 . . 68 (9 . . × − ( − . . . 1) 8 . 15 12 . . . 74 (7 . . 64 (6 . . × − ( − . . − ) 5 . 12 6 . 62 (29%) 6 . 62 (29%) 5 . 61 (9 . . 61 (9 . × − ( − . . . . 18 6 . 67 (29%) 6 . 63 (28%) 5 . 66 (9 . . 66 (9 . × − ( − . . . 1) 9 . 27 13 . . . 85 (6 . . 81 (5 . × − (0 . . − ) 4 . 23 5 . 46 (29%) 5 . 46 (29%) 4 . 64 (9 . . 64 (9 . . × − (0 . . . . 53 5 . 76 (27%) 5 . 71 (26%) 4 . 91 (8 . . 91 (8 . . × − (0 . . . 1) 19 . . . . . . . . . × − (0 . . − ) 5 . 12 6 . 62 (29%) 6 . 62 (29%) 5 . 61 (9 . . 61 (9 . . × − (0 . . . . 42 6 . 92 (27%) 6 . 88 (27%) 5 . 89 (8 . . 89 (8 . . × − (0 . . . 1) 21 . . . . . . . . . × − x = (100 , − T ( − . . − ) 3 . 72 4 . 96 (33%) 4 . 96 (33%) 4 . 13 (11%) 4 . 13 (11%) 4 . × − ( − . . . . 06 5 . 31 (31%) 5 . 26 (29%) 4 . 45 (9 . . 45 (9 . . × − ( − . . . 1) 19 . . . . . . . . × − ( − . . − ) 4 . 63 6 . 13 (33%) 6 . 13 (33%) 5 . 12 (10%) 5 . 12 (10%) 3 × − ( − . . . . 96 6 . 48 (30%) 6 . 42 (29%) 5 . 43 (9 . . 43 (9 . × − ( − . . . 1) 22 . . . . . . . . × − (0 . . − ) 3 . 72 4 . 96 (33%) 4 . 96 (33%) 4 . 13 (11%) 4 . 13 (11%) 4 . × − (0 . . . . 81 5 . 05 (32%) 5 . 01 (31%) 4 . 22 (10%) 4 . 22 (10%) 4 . × − (0 . . . 1) 9 . 93 19 . . . . . . . × − (0 . . − ) 4 . 62 6 . 13 (33%) 6 . 13 (33%) 5 . 12 (10%) 5 . 12 (10%) 2 . × − (0 . . . . 72 6 . 22 (32%) 6 . 18 (31%) 5 . 20 (10%) 5 . 20 (10%) 2 . × − (0 . . . 1) 11 . . . . . . . . × − Table 3: The cost of trading incurred for different values of x , ρ , ¯ η and λ , for the two asset case ( n = 2). The percentage values denote the increasecompared to the optimal solution. Smaller percentages are better. These results are based on 200 realizations of ξ to ξ . h table:ex2A i eferences h alfonsi2010optimal i A. Alfonsi, A. Schied, and A. Fruth. Optimal execution strategies in limit order books with generalshape functions. Quantitative Finance , 10(2):143–157, 2010. h alfonsi2009order i A. Alfonsi, A. Schied, and A. Slynko. Order book resilience, price manipulation, and the positive portfolioproblem. SIAM Journal on Financial Mathematics , 3(1):511–533, 2012. h almgren2003optimal i R. Almgren. Optimal execution with nonlinear impact functions and trading-enhanced risk. AppliedMathematical Finance , 10(1):1–18, 2003. h almgren2012optimal i R. Almgren. Optimal trading with stochastic liquidity and volatility. SIAM Journal on FinancialMathematics , 3(1):163–181, 2012. h almgren2001optimal i R. Almgren and N. Chriss. Optimal execution of portfolio transactions. Journal of Risk , 3:5–40, 2001. h AL2006 i R. Almgren and J. Lorenz. Bayesian adaptive trading with a daily cycle. Journal of Trading , 1(4):38–46,2006. h bertsimas1998optimal i D. Bertsimas and A. Lo. Optimal control of execution costs. Journal of Financial Markets , 1(1):1–50,1998. h bhatia2001linear i R. Bhatia. Linear algebra to quantum cohomology: the story of Alfred Horn’s inequalities. The AmericanMathematical Monthly , 108(4):289–318, 2001. h boogert2008gas i A. Boogert and C. De Jong. Gas storage valuation using a Monte Carlo method. The Journal ofDerivatives , 15(3):81–98, 2008. h easley1987price i D. Easley and M. O’hara. Price, trade size, and information in securities markets. Journal of FinancialEconomics , 19(1):69–90, 1987. h eisler2012price i Z. Eisler, J.-P. Bouchaud, and J. Kockelkoren. The price impact of order book events: market orders,limit orders and cancellations. Quantitative Finance , 12(9):1395–1419, 2012. h EF2007 i R. Engle and R. Ferstenberg. Execution risk. Journal of Portfolio Management , 33(2):34–44, 2007. h gatheral2010no i J. Gatheral. No-dynamic-arbitrage and market impact. Quantitative Finance , 10(7):749–759, 2010. h GSA2012 i J. Gatheral, A. Schied, and A. Slynko. Transient linear price impact and Fredholm integral equations. Mathematical Finance , 22:445–474, 2012. h HJ1985 i R. Horn and C. R. Johnson. Matrix Analysis . Cambridge University Press, 1985. h huberman2004price i G. Huberman and W. Stanzl. Price manipulation and quasi-arbitrage. Econometrica , 72(4):1247–1275,2004. h JK2008 i S. Joe and F. Kuo. Constructing Sobol’ sequences with better two-dimensional projections. SIAMJournal of Scientific Computing , 30(5):2635–2654, 2008. h karatzas1991brownian i I. Karatzas and S. Shreve. Brownian motion and stochastic calculus , volume 113. Springer Verlag, 1991. h konishiy2001optimal i H. Konishi. Optimal slice of a VWAP trade. Journal of Financial Markets , 5(2):197–221, 2002. h kyle1985continuous i A. Kyle. Continuous auctions and insider trading. Econometrica , 53(6):1315–1335, 1985. h longstaff2001valuing i F. Longstaff and E. Schwartz. Valuing American options by simulation: a simple least-squares approach. Review of Financial Studies , 14(1):113–147, 2001. h obizhaeva2005optimal i A. Obizhaeva and J. Wang. Optimal trading strategy and supply/demand dynamics. Journal of FinancialMarkets , 16(1):113–140, 2013. h potters2003more i M. Potters and J.-P. Bouchaud. More statistical properties of order books and price impact. Physica A:Statistical Mechanics and its Applications , 324(1):133–140, 2003. h schoneborn2011adaptive i T. Sch¨oneborn. Adaptive basket liquidation.