Dynamic Optimal Portfolios for Multiple Co-Integrated Assets
aa r X i v : . [ q -f i n . P M ] A ug Dynamic Optimal Portfolios for Multiple Co-Integrated Assets
T. N. L ˇi ∗ and A. Papanicolaou † August 7, 2019
Abstract
In this paper we construct and analyse a multi-asset model to be used for long-term sta-tistical arbitrage strategies. A key feature of the model is that all assets have co-integration ,which, if sustained, allows for long-term positive profits with low probability of losses. Optimalportfolios are found by solving a Hamilton-Jacobi-Bellman equation, to which we can intro-duce portfolio constraints such as market neutral or dollar neutral. Under specific conditions ofthe parameters, we can prove there is long-term stability for an optimal portfolio with stablegrowth rate. Historical prices of the S&P500 constituents can be tested for co-integration andour model calibrated for analysis, from which we find that co-integration strategies require aterminal investment horizon sufficiently far into the future in order for the optimal portfolios togain from co-integration. The data also demonstrates that statistical arbitrage portfolios willhave improved in-sample Sharpe ratios compared to multivariate Merton portfolios, and thatstatistical arbitrage portfolios are naturally immune to market fluctuations.
Keywords: co-integrated assets, market-neutral portfolios, statistical arbitrage, stochasticoptimal control, matrix Riccati equations.
AMS Subject Codes:
Statistical arbitrage strategies can be constructed by trading among pairs of assets having co-integrated prices. The essential idea is that a pair of co-integrated prices will have a difference thatis mean reverting. This mean-reverting difference is referred to as a spread . For an arbitrageur withthe financial means, a possible strategy is to long the cheaper asset, short the expensive asset, andthen to wait for the spread to converge, at which point the position can be closed for a profit. Thisis an example of statistical arbitrage because, while it may seem like a sure profit, there is no finitetime by when a spread will have almost-surely converged. Instead there is only a high probabilityof the spread converging before reaching a fixed, finite investment horizon.In this paper we construct a model for trading of multiple assets where each asset is co-integratedwith respect to a common benchmark index. Co-integration is found utilising the method of[9] where an asset’s log price is regressed against the log of the benchmark, and co-integrationconclusively identified if a unit-root hypothesis is rejected for the residual time series. As assumedin [19, 20], these residuals form a stationary Ornstein-Uhlenbeck (OU) process that henceforth is ∗ Department of Mathematics, New York University, 251 Mercer Street, New York, NY, 10012. [email protected] † Department of Finance and Risk Engineering, New York University, 1 MetroTech Center, Brooklyn NY, 11201. [email protected] W t , and the spreads Z it = a i + b i t + log (cid:0) S t (cid:1) + β i log (cid:0) S it (cid:1) for i = 1 , , · · · , d , where S t is the benchmark index, S it is the price of the i th tradedassets, and ( a i , b i , β i ) are the regression coefficients returned by the Engle-Granger test making Z it stationary. The optimal portfolio is the solution to a Hamilton-Jacobi-Bellman (HJB) partialdifferential equation (PDE), which in the case of power utility we are able to reduce to a system ofordinary differential equations (ODEs) that includes a matrix Riccati equation and a pair of linearequations. We perform long-term stability analysis of these ODEs, which gives us an indicationof the model’s soundness for financial application. In particular, finiteness of the ODEs for anyfinite-time investment horizon indicates that the model has an absence of arbitrage; if there wasan arbitrage then the ODEs would have a singularity. In addition, if the solution of the HJBconverges to an equilibrium as the investment horizon tends toward infinity, then there is a long-term statistical arbitrage portfolio that earns positive profits with probability close to one. Ouranalysis shows that how optimal strategies depend on the co-integration spreads, mean-reversionspeeds, risk aversion, and different portfolio constraints.We also perform empirical analysis based on historical financial market data. The purpose ofthese studies is to show that indeed, market data agree with our theoretical conditions for long-term stability. These studies also give us a sense of the effectiveness of implementing the optimalstrategies of these models. For our studies we select the subset of S&P500 constituents that areco-integrated with the S&P500 Index itself during the period of 2012-01-01 and 2019-05-15, andthen look at the long-term statistics of a portfolio of the 10 stocks whose spreads have the fastestrate of mean reversion. We find 106 co-integrated constituents, for which we then compute meanand covariance matrices of their returns processes, and then allocating among the top 10 fastestmean reverters we look at long-term Sharpe ratio for a few different optimal strategies. Thesestudies show (1) that statistical arbitrage strategies require a sufficiently-long investment horizonin order for the optimal portfolios to have a high probability of gaining from co-integration, and(2) that there is in-sample out-performance of co-integration-based portfolios in comparison totrading strategies based on constant-parameters with no consideration given to spreads (in otherwords, a multivariate Merton portfolio), and that statistical arbitrage portfolios have some naturalimmunity to market fluctuations such as the Chinese currency devaluation of August 2015. Westress that these data findings are in-sample tests, and merely serve to provide intuition and analysisfor statistical arbitrage opportunities; in-sample tests suffer from over-fitting and there are majorhurdles to implement out-of-sample tests because portfolio performance is significantly effected byparameter-estimation error. 2 .1 Overview of Related Work A test for co-integration of financial time series was designed in [9], namely, the Engle-Granger test.An application of [9] and some basic examples of co-integration-based trading rules were shown in[21], trading of co-integrated pairs alongside methods for filtering and parameter estimation tohandle latency was studied in [8], and an in-depth statistical analysis of the performance of pairstrading strategies was done in [13]. Principal component analysis of large number of assets co-integrated through common factors was the topic in [4], and is the basis for the model in this paper;empirical testing of pairs trading, including out-of-sample experiments with changing parameters,was completed in [12]. Analysis showing significance of short-term reversal and momentum factorson returns of pairs trading portfolios was finished in [5].The literature on stochastic control and optimisation for co-integration models also has somedepth. Stochastic optimal control for pairs trading with OU spreads was done in [19]. A stochasticcontrol for optimal trading of co-integrated pairs is posed and solved in [3, 20], and stochasticoptimal control for pairs trading with a local-volatility model was analysed by [18]. Optimaltrading of spreads with transaction costs and stop-loss criterion were solved and analysed in [17,16]. Related to the stability analysis presented in this paper, there was the long-term stabilityanalysis for matrix Riccati equations of multi-asset models done in [7], and the matrix Riccatiequations analysis for a single co-integrated pair done in [15]. There were also machine learningapproaches taken to statistical arbitrage, such as reinforcement learning and boosting applied toco-integrated constituents in the S&P500 were finished in [10]. An HJB equation for an optimalportfolio constrained to be 100% long was completed in [2] with application toward comparingactive and passive fund management.
In this paper, we presents and solve the stochastic control and optimisation problems and thenanalyse the solutions – both for the cases of unconstrained and constrained portfolios. Later onin the paper we perform some empirical studies on historical data. The organisation of the pa-per is as follows: Section 2 contains the definitions for the model along with analysis of the HJBequations, with Section 2.2 presenting the solution of the HJB for the unconstrained portfoliovia an exponential-affine ansatz, with Section 2.3 presenting the stability analysis, and then withSection 2.4 presenting the HJB and stability analysis for optimisation with market neutral anddollar-neutral constraints; the empirical analysis of historical data comes in Section 3 with pre-liminary data analysis and parameter estimation presented in Section 3.1 and analysis of portfolioperformance (e.g., Sharpe ratios) presented in Section 3.2; Section 4 is the conclusion.
This section introduces the stochastic model for multiple co-integrated assets, derives the HJBequations for optimal portfolios, and shows the stability analysis for the solutions.
Let (cid:16) Ω , F , ( F t ) t ≥ , P (cid:17) be a filtered probability space. Suppose S t is a benchmark for the financialmarket or an industrial sector with no uncertainty in expected returns. The stochastic differential3quation for its dynamic is d log S t = (cid:18) µ − σ (cid:19) dt + σ dB t , (2.1)where µ and σ > B t is an one dimensional Standard Brownian Motion on thefiltered probability space. Suppose S it , i = 1 , , · · · , d ∈ N , is the price of an individual firm in thefinancial market or the industrial sector, whose dynamic is modelled by the stochastic differentialequation d log S it = (cid:18) µ i − σ i δ i Z it (cid:19) dt + σ i dB it , (2.2)where µ i , δ i , and σ i > B it is an element of the vector B t = (cid:2) B t , B t , · · · , B dt (cid:3) ⊤ ∈ R d , which is a d -dimensional Standard Brownian Motions with correlation coefficient ρ ij ∈ ( − , i, j = 1 , , · · · , d such that dB it dB jt = ρ ij dt , and where (cid:2) B t , B it (cid:3) is a 2-dimensional StandardBrownian Motions with correlation coefficient ρ i ∈ ( − ,
1) such that dB it dB t = ρ i dt for i =1 , , · · · , d . Furthermore, the co-integrated process Z it that appears in equation (2.2) is spreadgiven by Z it = a i + b i t + log S t + β i log S it , (2.3)where typically β i < Z it is thedifferential of (2.3) dZ it = (cid:18) b i + µ − σ β i µ i − β i σ i β i δ i Z it (cid:19) dt + σ dB t + β i σ i dB it (2.4)= − β i δ i − b i + µ + β i µ i − σ − β i σ i β i δ i − Z it ! dt + σ dB t + β i σ i dB it . We assume that each element Z it of the vector Z t = (cid:2) Z t , Z t , · · · , Z dt (cid:3) ⊤ ∈ R d is a one dimensionalOU process, which holds true if − β i δ i >
0, in other words if δ i >
0. We also assume that there isa risk-free asset, such as a money market account, with interest rate r > Proposition 2.1.
Given correlation matrix ρ = [ ρ ij ] ∈ R d × d , i, j = 1 , , · · · , d , and correlationcolumn vector ~ ρ = (cid:2) ρ , ρ , · · · , ρ d (cid:3) ⊤ ∈ R d , in order for them to combine to form a ( d + 1) -dimensional correlation structure, they must satisfy the condition of ~ ρ ⊤ ρ − ~ ρ ≤ . Proof. Consider a linear representation of B t B t = ω ˜ B t + d X i =1 ω i B it , where ω ∈ R and ω i ∈ R , ˜ B t and B it are two independent Standard Brownian Motions. Thequadratic variation of B t is dt = dB t dB t = ω d ˜ B t + d X i =1 ω i dB it ! = (cid:16) ω + ω ⊤ ρω (cid:17) dt , here ω = [ ω , ω , . . . , ω d ] ⊤ ∈ R d . Thus, it must be that ≤ ω ⊤ ρω = 1 − ω ≤ ; and the crossvariation of B t and B it is ρ i dt = dB it dB t = dB it ω d ˜ B t + d X j =1 ω j dB jt = d X j =1 ω j ρ ij dt . Hence, ρ = ρω . Therefore, ~ ρ ⊤ ρ − ~ ρ = ( ρω ) ⊤ ρ − ( ρω )= ω ⊤ ρω = 1 − ω ≤ . This inequality describes the relationship between the correlations within the assets and the corre-lations between the benchmark and assets.
Example 2.1.
Suppose that ρ = I , where I ∈ R d × d is an identity matrix. Then it must be that k ~ ρ k ≤ to ensure ~ ρ ⊤ ρ − ~ ρ ≤ . Example 2.2.
Consider the d -dimensional case in which the correlation coefficients ~ ρ and ρ havethe following structure ρ = (1 − c ) I + c ⊤ ~ ρ = c , where = [1 , , · · · , ⊤ ∈ R d , c ∈ R , and c ∈ R with < c < . From the Sherman-Woodbury-Morrison formula, we have ρ − = 11 − c (cid:18) I − c d − c ⊤ (cid:19) . Therefore, ~ ρ ⊤ ρ − ~ ρ = c ⊤ ρ − = c d − c ≤ , which is the case if c ≤ (cid:0) − n (cid:1) c + n . The wealth process W t of an arbitrageur is constructed in the following way. We denote by π it the fractions of wealth invested in the i th risky asset at time t , then (cid:16) − P di =1 π it (cid:17) is the fractionsof wealth invested in the risk-free asset. Therefore, consider a wealth process that is among theclass of self-financing strategies and hence, its evolution is given by dW t = d X i =1 π it W t dS it S it + r − d X i =1 π it ! W t dt . (2.5)5lease note that the benchmark S t is not contained in the portfolio, as you can observe from theequation (2.5).Next, we define the value function u ( t, w, z ) for any ( t, w, z ) ∈ [0 , T ] × R + × R d where T < ∞ is a terminal time. The arbitrageur seeks a control π : [0 , T ] × Ω → R d where at each timethe portfolio allocation is π t = (cid:2) π t , π t , · · · , π dt (cid:3) ⊤ ∈ R d , to maximise the expectation of terminalutility, u ( t, w, z ) = sup π ∈A E (cid:2) U ( W T ) (cid:12)(cid:12) W t = w, Z t = z (cid:3) , (2.6)where w and z = [ z , · · · , z d ] ⊤ ∈ R d are the state variables, and π is selected from a set ofadmissible controls A defined as A = ( π : [0 , T ] × Ω → R d (cid:12)(cid:12)(cid:12) π t ∈ F t , and Z T k π t W t k dt < ∞ a.s. ) . (2.7)For a mathematical primer on the theory of stochastic control and optimisation, we refer the readerto chapter 4 of [11].We assume that the arbitrageur has a concave utility function U ( x ) over power type: U ( w ) = 1 γ w γ , (2.8)where γ < γ approaches to zero, it indicates that the arbitrageur is more risk loving, and if γ tends to −∞ , itrepresents that the arbitrageur is more risk averse. In order to have convenience in mathematical derivations for the rest sections of this paper, weintroduce some notations initially. For i, j = 1 , , · · · , d ∈ N , we denote the following vectors andmatrices: ~ µ = [ µ i ] ∈ R d , µ = [ µ i − r ] ∈ R d , (2.9) ~ σ = [ σ i ] ∈ R d , σ = diag ( ~ σ ) ∈ R d × d ,~ β = [ β i ] ∈ R d , β = diag (cid:16) ~ β (cid:17) ∈ R d × d ,~ δ = [ δ i ] ∈ R d , δ = diag (cid:16) ~ δ (cid:17) ∈ R d × d ,~ κ = [ − β i δ i ] ∈ R d , κ = diag ( ~ κ ) ∈ R d × d , b = [ b i ] ∈ R d , θ = " − b i + µ + β i µ i − (cid:0) σ + β i σ i (cid:1) β i δ i ∈ R d , Σ = [ σ i σ j ρ ij ] ∈ R d × d , Σ = (cid:2) σ σ i ρ i + σ i σ j β j ρ ij (cid:3) ∈ R d × d , Σ = (cid:2) σ + σ σ i β i ρ i + σ σ j β j ρ j + σ i β i σ j β j ρ ij (cid:3) ∈ R d × d . Please note that we assume that Σ is positive definite and invertible in this paper. We can alsoobserve that Σ and Σ are symmetric matrices, Σ , however, is not. Subsequently, we apply the6tandard stochastic control and optimisation techniques and expect the value function u defined in(2.6) to satisfy the following HJB equation: − u t − ( θ − z ) ⊤ κ ∇ z u −
12 tr (cid:0) Σ ∇ z u (cid:1) − rwu w (2.10) − sup π ∈ R d (cid:18) π ⊤ ( µ + δ z ) wu w + π ⊤ Σ ∇ z ( ∇ w u ) w + 12 π ⊤ Σ π w u ww (cid:19) = 0 ,u ( T, w, z ) = 1 γ w γ . The wealth variable w can be factored out of the solution by utilising the ansatz u ( t, w, z ) = 1 γ w γ g ( t, z ) , (2.11)hence, we have u t = 1 γ w γ g t , u ww = ( γ − w γ − g , ∇ z u = 1 γ w γ ∇ z g ,u w = w γ − g , ∇ z u = 1 γ w γ ∇ z g , ∇ z ( ∇ w u ) = w γ − ∇ z g . Therefore, the HJB equation (2.10) can be transformed into − g t − ( θ − z ) ⊤ κ ∇ z g −
12 tr (cid:0) Σ ∇ z g (cid:1) − rγg (2.12) − inf π ∈ R d (cid:18) π ⊤ ( µ + δ z ) γg + π ⊤ γ Σ ∇ z g + 12 π ⊤ Σ π gγ ( γ − (cid:19) = 0 ,g ( T, z ) = 1 . We compute the optimal control variable π ∗ by solving the optimisation problem that is containedinside the HJB equation (2.12) in terms of g and its partial derivatives π ∗ = 11 − γ Σ − ( µ + δ z ) + 11 − γ Σ − Σ ∇ z (log g ) . (2.13)Inserting the optimal π ∗ back into the HJB equation (2.12) results in the following nonlinear partialdifferential equation g t + ( θ − z ) ⊤ κ ∇ z g + 12 tr (cid:0) Σ ∇ z g (cid:1) + rγg (2.14)+ γg − γ ) ( µ + δ z + Σ ∇ z (log g )) ⊤ Σ − ( µ + δ z + Σ ∇ z (log g )) = 0 ,g ( T, z ) = 1 . We approach solving PDE (2.14) with an exponential-affine ansatz for g , g ( t, z ) = exp (cid:16) a ( t ) + b ⊤ ( t ) z + z ⊤ C ( t ) z (cid:17) , (2.15)where a ( t ) ∈ R is a scalar, b ( t ) = [ b i ( t )] ∈ R d is a column vector, and C ( t ) = [ c ij ( t )] ∈ R d × d , i, j = 1 , , · · · , d , is a d × d symmetric matrix. By utilising the ansatz (2.15), the nonlinearity inPDE (2.14) can be removed, and the equation can be transformed into a system of ODEs.7 roposition 2.2. The PDE (2.14) is solved utilising the exponential-affine an-satz of (2.15) , andthe functions a ( t ) ∈ R , b ( t ) ∈ R d , and C ( t ) ∈ R d × d satisfy the following system of ODEs: ∂a ( t ) ∂t = − b ⊤ ( t ) (cid:18) γ − γ ) Σ ⊤ Σ − Σ + 12 Σ (cid:19) b ( t ) (2.16) − γ − γ ) b ⊤ ( t ) Σ ⊤ Σ − µ − γ − γ ) µ ⊤ Σ − Σ b ( t ) − γ − γ ) µ ⊤ Σ − µ − rγ − θ ⊤ κb ( t ) − tr [ Σ C ( t )] ,a ( T ) = 0 ; ∂ b ( t ) ∂t = − C ( t ) (cid:18) γ − γ Σ ⊤ Σ − Σ + Σ (cid:19) b ( t ) (2.17) − (cid:18) γ − γ δ Σ − Σ − κ (cid:19) b ( t ) − C ( t ) (cid:18) γ − γ Σ ⊤ Σ − µ + κθ (cid:19) − γ − γ δ Σ − µ , b ( T ) = ; ∂ C ( t ) ∂t = − C ( t ) (cid:18) γ − γ Σ ⊤ Σ − Σ + 2 Σ (cid:19) C ( t ) (2.18) − C ( t ) (cid:18) γ − γ Σ ⊤ Σ − δ − κ (cid:19) − (cid:18) γ − γ Σ ⊤ Σ − δ − κ (cid:19) ⊤ C ( t ) − γ − γ ) δ Σ − δ , C ( T ) = . Proof. By inserting the exponential-affine ansatz (2.15) into PDE (2.14) , and grouping terms aseither quadratic in z , linear in z or constant in z , then equations (2.16) , (2.17) and (2.18) arerespectively obtained. We can observe that the three ODEs (2.16), (2.17) and (2.18) with respect to a ( t ), b ( t ) and C ( t ) are coupled. The coupling is recursive in the sense that the equation with respect to C ( t ) isautonomous, while we can solve for b ( t ) given ( C ( t )) u ≥ t , and solve for a ( t ) given ( b ( t )) u ≥ t and( C ( t )) u ≥ t . The type of the ODE that C ( t ) solves is a matrix Riccati equation. The gradient ofthe exponential-affine ansatz is ∇ z (log g ) = 2 C ( t ) z + b ( t ) so that the optimal π from equation(2.13) can be expressed as π ∗ = 11 − γ Σ − ( µ + δ z ) + 11 − γ Σ − Σ (2 C ( t ) z + b ( t )) . Verification that this is indeed an optimal strategy is straight-forward under the multiple-co-integrated model that we’ve constructed: 8 roposition 2.3 (Verification) . The optimal control variable π ∗ that is given by formula (2.13) utilising the exponential-affine ansatz that is proposed by (2.15) , where C ( t ) , b ( t ) and a ( t ) arerespectively the solutions of the matrix Riccati equation (2.18) , the ODEs (2.17) and (2.16) , belongsto the set of admissible controls A that is described by formula (2.7) and maximises the expectedutility function that is defined in formula (2.6) .Proof. The form of the model fits into the framework set forth in [6] and [7], and hence the sameargument for verification applies here. Stability analysis of the matrix Riccati equation (2.18) and the linear ODE (2.17) informs us thatwhether our solution to PDE (2.14) blows up or not. We extend the time domain for the ODEs(2.18) and (2.17) to ( −∞ , T ] for any finite T , and if the solution remains finite at all time then wehave a stable system from which we can draw intuition about long-term investment strategies.Our analysis of the matrix Riccati equation C ( t ) will use Theorem 2.1 from [23] to show thatthe solution of equation (2.18) exists, is bounded and unique for all t ≤ T . Let us rewrite the matrixRiccati equation (2.18) as ∂ C ( t ) ∂t = − A ⊤ u C ( t ) − C ( t ) A u − C ( t ) Q u C ( t ) − P u , (2.19) C ( T ) = , where Q u = 2 γ − γ Σ ⊤ Σ − Σ + 2 Σ , A u = γ − γ Σ ⊤ Σ − δ − κ , P u = γ − γ ) δ Σ − δ . Proposition 2.4.
For γ < , the coefficient matrices of the quadratic term and the zero order termof the matrix Riccati equation C ( t ) , namely, Q u and P u in equation (2.19) are symmetric positivedefinite and symmetric negative definite, respectively.Proof. From formula (2.9) , we have Σ = σ σ ~ ρ ⊤ + Σ β , (2.20) Σ = σ ⊤ + σ (cid:18) βσ ~ ρ ⊤ + (cid:16) ~ ρ ⊤ (cid:17) ⊤ βσ (cid:19) + β Σ β . Hence, the coefficient matrix Q u of the quadratic term of the matrix Riccati equation (2.19) hasthe following decomposition, Q u = 2 γ − γ (cid:16) σ σ ~ ρ ⊤ + Σ β (cid:17) ⊤ Σ − (cid:16) σ σ ~ ρ ⊤ + Σ β (cid:17) (2.21)+ 2 σ ⊤ + 2 σ (cid:18) βσ ~ ρ ⊤ + (cid:16) ~ ρ ⊤ (cid:17) ⊤ βσ (cid:19) = 21 − γ Σ − γσ − γ (cid:18) ⊤ − γ − γ ρ ⊤ ρ − ρ (cid:19) , here ρ = ~ ρ ⊤ ∈ R d × d . The matrix Σ is symmetric positive definite by construction becauseit is the diffusion matrix of the OU processes (2.4) , and it is straightforward to check that matrix ⊤ − γ − γ ρ ⊤ ρ − ρ is symmetric positive semidefinite for γ < . Hence, the third line of equation (2.21) is a symmetric positive definite matrix, and it follows that the coefficient matrix of thequadratic term of the matrix Riccati equation (2.19) is symmetric positive definite.Proving that P u is symmetric negative definite is uncomplicated. We can observe that matrix δ Σ − δ is symmetric positive definite. Consequently, for γ < , − P u is symmetric positive definite,in other words P u is symmetric negative definite. Given Proposition 2.4, the stability analysis from [23] applies directly, but first it will be usefulto define the following properties:
Definition 2.1 (Controllability) . Let A ∈ R n × n and B ∈ R n × m be constant matrices. The con-trollability matrix of ( A , B ) is the n × mn matrix Γ ( A , B ) = (cid:2) B , AB , · · · , A n − B (cid:3) . The pair ( A , B ) is controllable if the rank of Γ is n . If ( A , B ) is controllable, so is ( A − BM , B ) for every matrix M ∈ R m × n . Definition 2.2 (Observability) . Let A ∈ R n × n and E ∈ R p × n be constant matrices. The pair ( E , A ) is observable if the pair (cid:0) A ⊤ , E ⊤ (cid:1) is controllable. Definition 2.3 (Stabilizability) . Let A ∈ R n × n and B ∈ R n × m be constant matrices. The pair ( A , B ) is stabilizable if there exists a constant matrix M such that all the eigenvalues of A − BM have negative real parts. With these definitions we have the following proposition for the matrix Riccati equation (2.18).
Proposition 2.5.
For γ < , the coefficient matrix Q u of the quadratic term in the matrix Riccatiequation (2.19) is symmetric positive definite. Consequently, there are matrices B u , E u , and N u such that Q u = B u N − u B ⊤ u and − P u = E ⊤ u E u , with the pair ( A u , B u ) being stabilizable, andthe pair ( E u , A u ) being observable. Hence, there is an unique solution C ( t ) to matrix Riccatiequation (2.19) that is negative semidefinite and bounded on ( −∞ , T ] , and there exists a uniquelimit ¯ C = lim t →−∞ C ( t ) .Proof. We first perform a change of variable. Define ˜ C ( t ) = − C ( t ) , so the matrix Riccati equation (2.19) becomes ∂ ˜ C ( t ) ∂t + A ⊤ u ˜ C ( t ) + ˜ C ( t ) A u − ˜ C ( t ) Q u ˜ C ( t ) + ( − P u ) = 0 , (2.22)˜ C ( T ) = 0 . Proposition 2.4 has shown that matrix − P u ∈ R d × d is symmetric positive definite. So, all eigen-values λ P u of − P u are positive and there exists an orthonormal basis for R d of their associatedeigenvectors, in other words, there is an orthonormal matrix O u such that − P u = O u D u O ⊤ u ,where D u = diag ( λ P u ) ∈ R d × d is a diagonal matrix with positive entries on the diagonal. Hence,we can write − P u = E ⊤ u E u , where E u = (cid:0) O u √ D u (cid:1) ⊤ is a real square matrix. The matrix E u isinvertible, and so the controllability matrix Γ (cid:0) A ⊤ u , E ⊤ u (cid:1) ∈ R d × d , as defined by Definition 2.1, hasrank d . Consequently, the pair (cid:0) A ⊤ u , E ⊤ u (cid:1) is controllable and the pair ( E u , A u ) is observable as perDefinition 2.2. he symmetric positive definiteness of Q u is proven in Proposition 2.4 as well. Thus, the matrix Q u ∈ R d × d also has diagonal decomposition, Q u = B u N − u B ⊤ u , where B u is an orthogonal matrix,and N − u = diag ( λ Q u ) ∈ R d × d where λ Q u are the positive eigenvalues of Q u . The matrix B u isinvertible, and therefore we can find a constant matrix M u ∈ R d × d such that all eigenvalues of A u − B u M u have negative real parts, and therefore the pair ( A u , B u ) stabilizable.The above analysis of matrices Q u and − P u confirms that we can apply Theorem 2.1 from [23]to conclude that solution ˜ C ( t ) to the matrix Riccati equation (2.22) is unique, positive semidefinite,bounded on ( −∞ , T ] , and has unique limit as t tends toward −∞ . Example 2.3.
Suppose instead of utilising the power utility function (2.8) , we consider the expo-nential utility function U ( x ) = − exp ( − γx ) , where γ > . Then the ansatz function (2.11) becomes u ( t, w, z ) = − exp ( − γw ) g ( t, z ) . In the case, the HJB equation is g t + ( θ − z ) ⊤ κ ∇ z g + 12 tr (cid:0) Σ ∇ z g (cid:1) + rγg −
12 ( µ + δ z + Σ ∇ z (log g )) ⊤ Σ − ( µ + δ z + Σ ∇ z (log g )) = 0 ,g ( T, z ) = 1 . Therefore, the condition for long-term stability of the matrix Riccati equation with respect to C ( t ) is that the matrix Σ ⊤ Σ − Σ + Σ has to be symmetric positive definite. In other words exponentialutility function has the same matrix Riccati equation as the power utility function with γ → −∞ . Remark 1.
The stability analysis of Proposition 2.5 is sufficient for there to be no arbitrage inthe model proposed by equations (2.1) and (2.2) . If there were an arbitrage then it would alwaysbe optimal to take additional positions in the arbitrage portfolio, hence causing the value functionto reach a Nirvana (see [15]), but there will be singularity if there is stability of the matrix Riccatiequation with respect to C ( t ) . After solving the matrix Riccati equation (2.18), we then can study the stability of the solutionto the ODE with respect to b ( t ), in other words equation (2.17). We start the analysis by intro-ducing the following lemma (a theorem from [22]) in regard to the eigenvalues of matrices (generalbackground for this lemma can be found in Chapter 1 of [14]). Lemma 2.1 (Theorem from [22]) . Let M be a d × d matrix and define the field of values S ( M ) := (cid:8) v ⊤ Mv | v is a vector such that v ⊤ v = 1 (cid:9) , which contains the eigenvalues of M .(a) Suppose M and M are two d × d matrices. If λ is an eigenvalue of M + M , then λ ∈ S ( M ) + S ( M ) = { λ + λ | λ ∈ S ( M ) , λ ∈ S ( M ) } ;(b) Suppose M and M are two d × d matrices, / ∈ S ( M ) and M − exits. If λ is aneigenvalues of M − M , then λ ∈ S ( M ) / S ( M ) = { λ / λ | λ ∈ S ( M ) , λ ∈ S ( M ) } ;(c) Suppose M is an arbitrary d × d matrix, M is symmetric positive semidefinite matrix. If λ is an eigenvalues of M M , then λ ∈ S ( M ) S ( M ) = { λ λ | λ ∈ S ( M ) , λ ∈ S ( M ) } . Proposition 2.6.
Let R u ( t ) be the coefficient matrix of homogeneous part of equation (2.17) . For γ < , if ⊤ ~ ρ = P di ρ i ≥ , then there exists t ∗ > −∞ such that R u ( t ) has all positive eigenvaluesfor t < t ∗ , and therefore the solution of ODE (2.17) has a finite steady state. roof. By observing the ODE (2.17) and utilising the notations given by formulae (2.9) and (2.20) ,we have the coefficient matrix for the ODE, R u ( t ) = − C ( t ) (cid:18) γ − γ Σ ⊤ Σ − Σ + Σ (cid:19) − (cid:18) γ − γ δ Σ − Σ − κ (cid:19) (2.23)= − C ( t ) Q u + 11 − γ κ − γσ − γ δσ − ρ − ρ . Let ¯ C = lim t →−∞ C ( t ) , then the limit of equation (2.23) is R u = − ¯ CQ u + 11 − γ κ − γσ − γ (cid:0) ρσδ − (cid:1) − ρ . (2.24) We can observe that matrices κ , δ , σ , and ρ are symmetric positive definite. Proposition 2.4 provesthat matrix Q u is symmetric positive definite, and proposition 2.5 implies that − ¯ C is symmetricpositive semidefinite. Hence, both κ and Q u are symmetric positive definite and so all the elementsof sets S ( κ ) and S ( Q u ) are positive, and − ¯ C is symmetric positive semidefinite and so all theelements of set S (cid:0) − ¯ C (cid:1) are non-negative. By the symmetric positive definiteness of those matrices,their product ρσδ − is a positive definite matrix as well. Also, because ρ = ~ ρ ⊤ , so rank (cid:0) ~ ρ ⊤ (cid:1) =1 , and ρ ~ ρ = ~ ρ ⊤ ~ ρ = (cid:16) ⊤ ~ ρ (cid:17) ~ ρ , it follows that ⊤ ~ ρ is an eigenvalue of ρ . Therefore, if ⊤ ~ ρ = P di ρ i ≥ , then ρ has non-negative eigenvalues.We then apply Lemma 2.1 on the matrix R u that is given by equation (2.24) . By item (c), we seethat matrix − ¯ CQ u has non-negative eigenvalues. By item (b), we see that matrix − γσ − γ (cid:0) ρσδ − (cid:1) − ρ has non-negative eigenvalues, and then by item (a) we see that − γ κ − γσ − γ (cid:0) ρσδ − (cid:1) − ρ has pos-itive eigenvalues. Hence, by item (a), we have that matrix R u has positive eigenvalues. Therefore,there exists a t ∗ > −∞ such that R u ( t ) has positive eigenvalues for all t < t ∗ , and the solution b ( t ) of the equation (2.17) has a finite steady state. Proposition 2.7.
The long-term growth rate of the certainty equivalent is proportional to thesolution a ( t ) of the ODE (2.16) .Proof. Given the utility function (2.8) , the value function (2.11) , and the expo- nential-affine ansatzfor g that is defined by formula (2.15) , the certainty equivalent is defined by U − ( u ( t, w, z )) = U − (cid:18) γ w γ g ( t, z ) (cid:19) = w exp (cid:18) γ (cid:16) a ( t ) + b ⊤ ( t ) z + z ⊤ C ( t ) z (cid:17)(cid:19) . Hence, under the optimal control variable π ∗ , the long-term growth rate is ln (cid:0) U − ( u ( t, w, z )) (cid:1) T − t = 1 T − t log (cid:18) w exp (cid:18) γ (cid:16) a ( t ) + b ⊤ ( t ) z + z ⊤ C ( t ) z (cid:17)(cid:19)(cid:19) = 1 γ ( T − t ) (cid:16) γ log ( w ) + b ⊤ ( t ) z + z ⊤ C z + a ( t ) (cid:17) . rom the analysis in Propositions 2.5 and 2.6, both solutions C ( t ) and b ( t ) to the ODEs (2.18) and (2.17) have finite limits as t tends to −∞ , and so it follows that a ( t ) is asymptotically linearas t tends to −∞ , and therefore, lim T →∞ T − t log (cid:0) U − ( u ( t, w, z )) (cid:1) (2.25)= lim t →−∞ γ ( T − t ) (cid:18) γ log ( w ) + b ⊤ ( t ) z + z ⊤ C ( t ) z + a ( t ) (cid:19) = lim t →−∞ γ ( T − t ) a ( t ) . Furthermore, denote by L ( t ) the right hand side of the ODE with respect to a ( t ) , in other wordsequation (2.16) . As T tends toward infinity we have lim t →−∞ L ( t ) = − ¯ b ⊤ (cid:18) γ − γ ) Σ ⊤ Σ − Σ + 12 Σ (cid:19) ¯ b − θ ⊤ κ ¯ b − tr (cid:0) Σ ¯ C (cid:1) − γ − γ ) ¯ b ⊤ Σ ⊤ Σ − µ − γ − γ ) µ ⊤ Σ − Σ ¯ b − γ − γ ) µ ⊤ Σ − µ − rγ = ¯ L , where ¯ b = lim t →−∞ b ( t ) and ¯ C = lim t →−∞ C ( t ) . Hence, as t tends toward negative infinity, equation (2.16) relaxes and we have lim t →−∞ T − t Z Tt ∂a ( s ) ∂s ds = ¯ L , and therefore the limit in equation (2.25) is lim T →∞ T − t log (cid:0) U − ( u ( t, w, z )) (cid:1) = 1 γ ¯ L , and the long-term growth rate is a constant.
It is common practice to seek statistical arbitrage strategies that have market neutrality. Market-neutrality generally means that the returns of a portfolio are impacted only by the idiosyncraticreturns of the assets contained in the portfolio, and are uncorrelated with the returns of a benchmarkor market factors. Hence, under the condition of market neutrality, if we can diversify with alarge number of co-integrated assets, then there is a very high probability that the portfolio canmaintain steady growth and low volatility. Market neutrality was discussed at length in the principalcomponent analysis of [4]. We will consider a portfolio to be market neutral if dW t W t dS t S t = 0. This13appens as follows, dW t W t dS t S t = d X i =1 π it dS it S it dS t S t = d X i =1 π it d log( S it ) (cid:0) dZ it − β i log S it (cid:1) = d X i =1 π it (cid:0) Σ ii − β Σ ii (cid:1) dt = σ d X i =1 π it σ i ρ i dt = 0 . where Σ ii and Σ ii are the i th elements of the diagonals of matrices Σ and Σ respectively.In matrix/vector notation, market neutrality is π ⊤ s = 0, where s = σ [ σ ρ , σ ρ , · · · , σ d ρ d ] ⊤ ∈ R d . If π ⊤ s = 1, then the portfolio has a unit-loading on the benchmark return. Dollar neutralityis π ⊤ = 0, where = [1 , , · · · , ⊤ ∈ R d , in other words the dollar amount invested in therisky assets sums to zero. If π ⊤ = 1 then the dollar amount invested in the riskless asset is zeroso that the portfolio is 100% long. Consequently, we obtain the following four different equalityconstraints: π ⊤ s = 0 : Constraint for market neutral , (2.26) π ⊤ s = 1 : Constraint for unit loading on market portfolio’s returns , π ⊤ = 0 : Constraint for dollar neutral , π ⊤ = 1 : Constraint for 100% long . With the constraints of (2.26) in mind, we now reformulate the optimal portfolio that is studiedin Section 2.1 and Section 2.2. We shall start from a general situation, in other words, we wantto integrate an equality constraint of the control variables π s + π s + · · · + π d s d = π ⊤ s = s π into the stochastic control and optimisation problem, which can include all the four different casesthat are given by formula (2.26), where s = [ s , s , · · · , s d ] ∈ R d is column vector and s π ∈ R is aconstraint coefficient. Therefore, the HJB equation (2.10) becomes − u t − ( θ − z ) ⊤ κ ∇ z u −
12 tr (cid:0) Σ ∇ z u (cid:1) − rwu w (2.27) − sup π ∈ R d π ⊤ s = s π (cid:18) π ⊤ ( µ + δ z ) wu w + π ⊤ Σ ∇ z ( ∇ w u ) w + 12 π ⊤ Σ π w u ww (cid:19) = 0 .u ( T, w, z ) = 1 γ w γ . Hence the transformed HJB equation (2.12) becomes − g t − ( θ − z ) ⊤ κ ∇ z g −
12 tr (cid:0) Σ ∇ z g (cid:1) − rγg (2.28) − inf π ∈ R d π ⊤ s = s π (cid:18) π ⊤ ( µ + δ z ) γg + π ⊤ γ Σ ∇ z g + 12 π ⊤ Σ π gγ ( γ − (cid:19) = 0 ,g ( T, z ) = 1 . λ denote a Lagrange multiplier, we define the Lagrangian function for the constrainedcontrol model L ( π , λ ) = π ⊤ ( µ + δ z ) γg + π ⊤ γ Σ ∇ z g + 12 π ⊤ Σ π gγ ( γ − − λ (cid:16) π ⊤ s − s π (cid:17) . Therefore, by the first order condition ∇ π L = , we can get the optimal values for the controlvariables π ∗ = 1 gγ (1 − γ ) Σ − ( gγ µ + gγ δ z + γ Σ ∇ z g − λ s ) (2.29)= 11 − γ Σ − (cid:18) µ + δ z + Σ ∇ z (log g ) − λgγ s (cid:19) . We then solve it for the Lagrange multiplier λ to get its optimal value with respect to the conditionof constraint s π = ( π ∗ ) ⊤ s = s ⊤ π ∗ λ ∗ = 1 s ⊤ Σ − s (cid:16) γ s ⊤ Σ − ( g µ + g δ z + Σ ∇ z g ) − gγ (1 − γ ) s π (cid:17) . (2.30)Inserting λ ∗ that is given by formula (2.30) back into equation (2.29), we can get the optimalcontrol variables π ∗ with respect to the constraint coefficient s π π ∗ = s π d c Σ − s + 11 − γ (cid:18) Σ − − d c Σ c (cid:19) ( µ + δ z ) (2.31)+ 11 − γ (cid:18) Σ − − d c Σ c (cid:19) Σ ∇ z (ln g ) , where d c = s ⊤ Σ − s ∈ R is a scalar, Σ c = Σ − ss ⊤ Σ − ∈ R d × d is a symmetric matrix. We theninsert the optimal control variable π ∗ that is given by formula (2.31) back into the HJB equation(2.28) and get the following non-linear partial differential equation g t + ( θ − z ) ⊤ κ ∇ z g + 12 tr (cid:0) Σ ∇ z g (cid:1) + rγg (2.32)+ gγ ( γ − s π d c + γgs π d c s ⊤ Σ − ( µ + δ z + Σ ∇ z (ln g ))+ γg − γ ) ( µ + δ z + Σ ∇ z (ln g )) ⊤ (cid:18) Σ − − d c Σ c (cid:19) × ( µ + δ z + Σ ∇ z (ln g )) = 0 g ( T, z ) = 1 . Similar to the unconstrained stochastic control problem of Section 2.2, we approach removingthe nonlinearity in PDE (2.32) by utilising the exponential-affine ansatz of (2.15), from which theoptimal portfolio can be expressed as π ∗ = s π d c Σ − s + 11 − γ (cid:18) Σ − − d c Σ c (cid:19) ( µ + δ z )+ 11 − γ (cid:18) Σ − − d c Σ c (cid:19) Σ (2 C ( t ) z + b ( t )) , where a ( t ), b ( t ) and C ( t ) are the solutions to the following system of ODEs:15 roposition 2.8. The PDE (2.32) is solved with the exponential-affine ansatz of (2.15) , and thefunctions a ( t ) ∈ R , b ( t ) ∈ R d , and C ( t ) ∈ R d × d satisfy the following system of ODEs: ∂a ( t ) ∂t = − b ⊤ ( t ) (cid:18) γ − γ ) Σ ⊤ (cid:18) Σ − − d c Σ c (cid:19) Σ + 12 Σ (cid:19) b ( t ) (2.33) − γ − γ ) µ ⊤ (cid:18) Σ − − d c Σ c (cid:19) µ − tr [ Σ C ( t )] − γ − γ ) b ⊤ ( t ) Σ ⊤ (cid:18) Σ − − d c Σ c (cid:19) µ − γ ( γ − s π d c − γ − γ ) µ ⊤ (cid:18) Σ − − d c Σ c (cid:19) Σ b ( t ) − rγ − (cid:18) θ ⊤ κ + γs π d c s ⊤ Σ − Σ (cid:19) b ( t ) − γs π d c s ⊤ Σ − µ ,a ( T ) = 0 ; ∂ b ( t ) ∂t = − C ( t ) (cid:18) γ − γ Σ ⊤ (cid:18) Σ − − d c Σ c (cid:19) Σ + Σ (cid:19) b ( t ) (2.34) − (cid:18) γ − γ δ (cid:18) Σ − − d c Σ c (cid:19) Σ − κ (cid:19) b ( t ) − C ( t ) (cid:18) γ − γ Σ ⊤ (cid:18) Σ − − d c Σ c (cid:19) µ + κθ (cid:19) − γ − γ δ (cid:18) Σ − − d c Σ c (cid:19) µ − γs π d c h C ( t ) Σ ⊤ + δ i Σ − s , b ( T ) = ; ∂ C ( t ) ∂t = − C ( t ) (cid:18) γ − γ Σ ⊤ (cid:18) Σ − − d c Σ c (cid:19) Σ + 2 Σ (cid:19) C ( t ) (2.35) − C ( t ) (cid:18) γ − γ Σ ⊤ (cid:18) Σ − − d c Σ c (cid:19) δ − κ (cid:19) − (cid:18) γ − γ Σ ⊤ (cid:18) Σ − − d c Σ c (cid:19) δ − κ (cid:19) ⊤ C ( t ) − γ − γ ) δ (cid:18) Σ − − d c Σ c (cid:19) δ , C ( T ) = . Proof. The proof is the same as that for Proposition 2.2.
Proposition 2.9 (Verification) . The optimal control variable π ∗ that is given by formula (2.13) utilising the exponential-affine ansatz that is proposed by formula (2.15) , where C ( t ) , b ( t ) and a ( t ) are respectively the solutions of matrix Riccati equation (2.35) , ODEs (2.34) and (2.33) , belongsto the set of admissible controls A that is described by formula (2.7) and maximises the expectedutility function that is defined in formula (2.6) .Proof. Similar to Proposition 2.3 for the unconstrained stochastic control problem, the verificationmethod of [6] and [7] applies. C ( t ) ofthe constrained stochastic control problem, which is expressed by equation (2.35). We rewrite thismatrix Riccati equation in the following way ∂ C ( t ) ∂t = − C ( t ) Q c C ( t ) − C ( t ) A c − A ⊤ c C ( t ) − P c , (2.36) C ( T ) = , where Q c = 2 γ − γ Σ ⊤ (cid:18) Σ − − d c Σ c (cid:19) Σ + 2 Σ , A c = γ − γ Σ ⊤ (cid:18) Σ − − d c Σ c (cid:19) δ − κ , P c = γ − γ ) δ (cid:18) Σ − − d c Σ c (cid:19) δ . Similar to the unconstrained stochastic control problem, we prove that the coefficient matrix ofits quadratic term is positive definite.
Proposition 2.10.
For γ < , the quadratic term Q c of the matrix Riccati equation (2.36) issymmetric positive definite and the matrix P c is symmetric negative semidefinite.Proof. Proposition 2.4 has shown that Q u is positive definite. By utilising formula (2.20) , Q c hasthe following decomposition Q c = 2 Σ + 2 γ − γ Σ ⊤ (cid:18) Σ − − d c Σ c (cid:19) Σ = 2 Σ + 2 γ − γ Σ ⊤ (cid:18) Σ − − d c Σ − ss ⊤ Σ − (cid:19) Σ = Q u − γd c (1 − γ ) Σ ⊤ Σ − ss ⊤ Σ − Σ . Because Σ is symmetric positive definite, so its inverse Σ − is symmetric positive definite as well.Hence d c = s ⊤ Σ − s > . Therefore, for γ < , Q c is positive definite.In order to prove P c is symmetric negative semidefinite, we first examine the symmetric ma-trix Σ − − d c Σ c . Utilising the Cauchy-Schwarz inequality ( x ⊤ Σ − s ) ≤ ( x ⊤ Σ − x )( s ⊤ Σ − s ) , weobserve that for any x ∈ R d , x ⊤ (cid:18) Σ − − d c Σ c (cid:19) x = x ⊤ Σ − x − x ⊤ Σ − ss ⊤ Σ − xs ⊤ Σ − s ≥ , where the equality holds if and only if x = s . Hence, Σ − − d c Σ c is positive semidefinite. Conse-quently, δ (cid:16) Σ − − d c Σ c (cid:17) δ is a symmetric positive semidefinite matrix with δ − s being the singlenull vector. Therefore, − P c is symmetric positive semidefinite with null-space spanned by δ − s . Similar to Proposition 2.5, because Q c is symmetric positive definite, − P c is symmetric positivesemidefinite, theorem 2.1 in [23] applies directly, and the solution to matrix Riccati equation (2.35)exists, is bounded and is unique. 17 roposition 2.11. For γ < and κ I , where I ∈ R d × d is the identity matrix, the coefficientmatrix Q c of the quadratic term in the matrix Riccati equation (2.36) of the constrained stochasticcontrol problem is symmetric positive definite and − P c is symmetric positive semidefinite. Con-sequently there are matrices B c , E c , and N c such that Q c = B c N − c B ⊤ c and − P c = E ⊤ c E c , withthe pair ( A c , B c ) being stabilizable, and the pair ( E c , A c ) being observable. Hence, there is anunique solution C ( t ) to matrix Riccati equation (2.36) that is negative semidefinite and boundedon ( −∞ , T ] , and there exists a unique limit ¯ C = lim t →−∞ C ( t ) .Proof. The symmetric positive semidefiniteness of the matrix − P c comes from Proposition 2.10,for which the proof shows that Σ − − d c (cid:0) Σ − s (cid:1) (cid:0) Σ − s (cid:1) ⊤ is symmetric positive semidefinite. Thismatrix is diagonalizable, Σ − − d c (cid:0) Σ − s (cid:1) (cid:0) Σ − s (cid:1) ⊤ = O c D c O ⊤ c , where O c ∈ R d × d is orthonormaland D c is a diagonal matrix with the non-negative eigenvalues of Σ − − d c (cid:0) Σ − s (cid:1) (cid:0) Σ − s (cid:1) ⊤ on itsdiagonal. Thus − P c = E ⊤ c E c , where E c = q − γ − γ ) √ D c O ⊤ c δ . As shown in the proof of Proposition2.10, there is an one dimensional null space of − P c , spanned by δ − s , and hence the rank of matrix E c is d − also with null space spanned by δ − s . However, the rank of the controllability matrix Γ (cid:0) A ⊤ c , E ⊤ c (cid:1) ∈ R d × d given by Definition 2.1 is d if k I . Indeed , recall the matrix A c definedfor the matrix Riccati equation (2.36) , A ⊤ c = γ − γ δ (cid:16) Σ − − d c Σ c (cid:17) ⊤ Σ − κ = 2 E ⊤ c E c δ − Σ − κ ,which when multiplied on the left-hand side by δ − s and on the right by E ⊤ c , (cid:16) δ − s (cid:17) ⊤ A ⊤ c E ⊤ c = (cid:16) δ − s (cid:17) ⊤ (cid:16) E ⊤ c E c δ − Σ − κ (cid:17) E ⊤ c = 2 (cid:16) E c δ − s (cid:17) ⊤ E c δ − Σ E ⊤ c − (cid:16) δ − s (cid:17) ⊤ κ E ⊤ c = − (cid:16) E c κδ − s (cid:17) ⊤ = 0 , if κ I . Therefore, if κ I , where I is the identity matrix, then Γ (cid:0) A ⊤ c , E ⊤ c (cid:1) ∈ R d × d has full rank. Thus,the pair (cid:0) A ⊤ u , E ⊤ u (cid:1) is controllable and the pair ( E u , A u ) is observable as per Definition 2.2.The symmetric positive definiteness of matrix Q c is proven in Proposition 2.10 as well. Similarto Proposition 2.5, Q c = B c N − c B ⊤ c , where B c is an orthogonal matrix, λ Q c are the eigenvalues of Q c , and N − c = diag ( λ Q c ) ∈ R d × d . The matrix B c is invertible, therefore we can find a constantmatrix M c ∈ R d × d such that all eigenvalues of matrix A c − B c M c have negative real parts, andtherefore the pair ( A c , B c ) is stabilizable.Finally, we let ˜ C ( t ) = − C ( t ) , so that the matrix Riccati equation (2.36) becomes ∂ ˜ C ( t ) ∂t + A ⊤ u ˜ C ( t ) + ˜ C ( t ) A u − ˜ C ( t ) B c N − c B ⊤ c ˜ C ( t ) + E ⊤ c E c = 0 , (2.37)˜ C ( T ) = 0 , and the above analysis of matrices Q c and − P c confirms that we can apply Theorem 2.1 from [23]again, to conclude that solution ˜ C ( t ) to equation (2.37) is unique, positive semidefinite, boundedon ( −∞ , T ] , and has unique limit as t tends toward −∞ . Next we analyse the behaviour of the solution for the ODE with respect to b ( t ) of the constrainedstochastic control problem that is described by equation (2.34).18 roposition 2.12. Let R c ( t ) be the coefficient matrix for the homogeneous part of equation (2.34) .For γ < , if ⊤ ~ ρ ≥ then there exists t ∗ > −∞ such that R c ( t ) has all positive eigenvalues for t < t ∗ . Therefore, the solution of the equation (2.34) has a finite steady state.Proof. By observing the equation (2.34) and following the notations that are denoted by formula (2.9) and formula (2.20) , we have R c ( t ) = − C ( t ) (cid:18) γ − γ Σ ⊤ (cid:18) Σ − − d c Σ c (cid:19) Σ + Σ (cid:19) − (cid:18) γ − γ δ (cid:18) Σ − − d c Σ c (cid:19) Σ − κ (cid:19) = − C ( t ) Q c + κ − γ − γ δ Σ − (cid:16) σ σ ~ ρ ⊤ + Σ β (cid:17) + γ − γ δ d c Σ c (cid:16) σ σ ~ ρ ⊤ + Σ β (cid:17) = − C ( t ) Q c + 11 − γ κ + 1 d c (1 − γ ) δ Σ − ss ⊤ γ β − γσ − γ δ (cid:18) Σ − − d c Σ − ss ⊤ Σ − (cid:19) σρ . Because of γ < and the assumption for the co-integrated vector, in other words β i < , the matrix γ β is positive definite. From Proposition 2.10, we know that matrix Q c is positive definite, fromProposition 2.11, we have that matrix − C ( t ) is positive semidefinite for t < t ∗ with t ∗ > −∞ ,and by the assumption of the model, we have that matrices κ and δ are positive definite. Notealso that Σ − − d c Σ − ss ⊤ Σ − is symmetric positive semidefinite. Hence, with the attribute for theeigenvalue of matrix ρ that is proven in Proposition 2.6, if ⊤ ~ ρ ≥ then by Lemma 2.1, R c ( t ) has all positive eigenvalues for t < t ∗ with t ∗ > −∞ . Therefore, the solution of the equation (2.34) is stable. Remark 2.
The long-term growth rate of the certainty equivalent of the constrained stochastic con-trol problem is proportional to the solution a ( t ) of equation (2.33) . The proof of such a propositionis similar to the proof of Proposition 2.7. This section describes the data, parameter estimation and numerical methods of estimating pa-rameters for solving the optimal controls. The results demonstrate how the formulae in Section 2apply to real-life finance.
We utilise the Yahoo Finance as our data source. The data set that we utilise is the adjusted dailyclose stock prices of the S&P 500 constituents from 2012-01-03 to 2019-05-15, and also includesthe SPY ETF among these traded assets. The non-traded benchmark S is the S&P500 Index.We assume that the interest rate r is 0.02, and in every calendar year, there are 252 trading days.The three parameters for the OU process Z it that is defined by formula (2.3) are ( a i , b i , β i ) for i = 1 , , · · · , d = 501, and are the outputs of the Engle-Granger co-integration test. By includingthe parameter a i , the process Z it has a stationary mean that is not significantly different from19ero. After finding all the co-integrated pairs, we select 10 out of the 106 companies that are thefastest mean reverters, in other words such that their values of − β i δ i are the largest ones. In orderto estimate the parameter δ i for the co-integrated process that is described by equation (2.4), weestimate an order-one auto-regressive coefficient for the discrete time series of the Z it process vialeast squares estimation, and then take logarithms, divide by ∆ t = 1 / − β i toobtain κ i . Table 3.1 shows the 10 fastest mean reverters alongside their Sharpe ratios and estimatedrates of mean reversion.For the 10 selected assets for trading, it remains to calculate their parameters for the returnsmodel. The expected rate of returns µ i are estimated with the a sample mean/median, but tendto be higher than prior views would suggest, so we place a multiplier of 0 . µ i = 0 . × mean[∆ S it /S it ]is such that 1% ≤ ˆ µ i ≤ Σ = [ σ i σ j ρ ij ] ∈ R × is estimated with astandard method of moments. Figure 3.1 shows the estimated Z it processes. Ticker Sharpe Ratio Mean-Reversion Rate
HST 0.338003 7.507944URI 0.651762 7.125851WDC 0.346716 7.624259APH 1.011030 8.724717MU 0.713108 9.123977AMP 0.698346 9.321913BWA 0.170926 9.305448TEL 0.780984 7.016104TIF 0.329176 7.433595SPY 0.915641 43.844862Table 3.1: Sharpe ratios and mean-reversion rates for the 10 fastest mean-reverting stocks. TheSharpe ratios are raw time-series estimates.
Process of Z t HSTURI WDCAPH MUAMP BWATEL TIFSPY
Figure 3.1: Co-integrated process Z t for the 10 fastest mean reverters.After estimating all parameters in the stochastic system, we then solve the matrix Riccatiequation for C ( t ) by applying Radon’s Lemma (see Chapter 3 of [1]). The method is a forward-time scheme, so we first perform the change of variable τ = T − t , and then the matrix Riccati20quation (2.22) for the unconstrained stochastic control problem becomes ∂ ˜ C ( τ ) ∂τ = A ⊤ u ˜ C ( τ ) + ˜ C ( τ ) A u − ˜ C ( τ ) Q u ˜ C ( τ ) − P u , ˜ C ( τ ) = 0 . Then the solution is given by ˜ C ( τ ) = ˜ C ( τ ) ˜ C − ( τ ) , where ˜ C ( τ ) is a solution of the initial value problem ∂ ˜ C ( τ ) ∂τ = (cid:16) − A u + Q u ˜ C ( τ ) (cid:17) ˜ C ( τ ) , ˜ C ( τ ) = ˜ C (0) , and h ˜ C ( τ ) , ˜ C ( τ ) i ⊤ is a solution of the associated linear system of ODEs ddt (cid:20) ˜ C ( τ )˜ C ( τ ) (cid:21) = (cid:20) − A u Q u − P u A ⊤ u (cid:21) (cid:20) ˜ C ( τ )˜ C ( τ ) (cid:21) , where ˜ C (0) is the initial condition, and matrices A u , P u , and Q u are given by formula (2.19). Tosolve the matrix Riccati equation (2.35) for the constrained stochastic control problem, we applythe same numerical method. Then we solve the ODEs with respect to b ( t ) and a ( t ) by utilisingthe fourth-order Runge-Kutta method.For γ = −
1, the numerical solutions are illustrated in Figure 3.2, Figure 3.3, and Figure 3.4(note that the time change-of-variable has been reversed so that the system is plotted backwardin time as the original problem had been posed). As we can observe clearly from these figures,the solutions C ( t ) and b ( t ) to matrix Riccati equation (2.18) and ODE (2.17) are stable, as thesufficient condition stated in Propositions 2.1, 2.4, 2.5, 2.6, 2.10, 2.11, and 2.12 are satisfied so longas γ < C ( t ), for unconstrainedstochastic control problem, utilising parameters estimated from the historical data. Notice thesystem starts to resemble its steady state if the trading time if there are 2 or more trading years. Ifthe trading time is less, then there is significant probability that the spreads may not converge tozero by the terminal time. Hence, the statistical arbitrage strategy will have a noticeable marginalimprovement if extra time is allotted. 21 O D E b ( t ) Asymptotic Behaviour Of Ordinary Differential Equation b(t) Of Model Without Constraint, = -1
Figure 3.3: The numerical solution of linear ODE (2.17) for b ( t ), for the unconstrained stochasticcontrol problem, utilising parameters estimated from the historical data. Notice the system startsto resemble its steady state if the trading time if there are 2 or more trading years. If the tradingtime is less, then there is significant probability that the spreads may not converge to zero by theterminal time. Hence, the statistical arbitrage strategy will have a noticeable marginal improvementif extra time is allotted. O D E a ( t ) Asymptotic Behaviour Of Ordinary Differential Equation a(t) Of Model Without Constraint, = -1
Figure 3.4: The numerical solution of linear ODE (2.18) for a ( t ), for the unconstrained stochasticcontrol problem, utilising parameters estimated from the historical data.Figures 3.2, 3.3, and 3.4 give us some idea of the required time for a statistical arbitrage strategyto be effective. From the plots showing the solutions of the ODEs, we see that the co-integratedsystem’s value function will need between one and two years before the marginal gain of extratrading time becomes log linear. If the allotted trading time for a statistical arbitrage is too short,then there is significant probability that spreads will not converge, and hence, for short time therewill noticeable marginal improvement if more time is allowed. After we solve the system of ODEs for the unconstrained and constrained stochastic control prob-lems, we then calculate historical wealth processes given by equation (2.5), and compare Sharpe22atios for different γ and different constraints. We work in the setting of very large terminal time T and consider a steady-state portfolio to demonstrate the results. In other words, we work withthe limiting vector π ∗ that is calculated utilising ¯ C = lim t →−∞ C ( t ) and ¯ b = lim t →−∞ b ( t ), e.g.,for the unconstrained stochastic control problem we have π ∗ ( z ) = 11 − γ Σ − (cid:0) µ + Σ ¯ b + ( δ + 2 Σ ¯ C ) z (cid:1) . Table 3.2 presents the annualised statistics for the optimal portfolio under the multivariate Mertonmodel with no co-integration (i.e., π merton = − γ Σ − µ ), Table 3.3 presents annualised statisticsfor the unconstrained portfolios under the model with co-integration, and Table 3.4 presents an-nualised statistics for the constrained portfolios under the model with co-integration. From thesetables we can see that as the risk aversion coefficient γ becomes more negative, in other words,as the arbitrageur becomes more risk averse, then there is a decrease in the expectation of excessreturn, the volatility of excess return, and the Sharpe ratio. Overall, the results of these numericalexperiments demonstrate that the optimal portfolio with co-integration performs significantly bet-ter than a constant-coefficient Merton portfolio. We can also observe that for the same values of γ ,the unconstrained portfolios in Table 3.3 have higher expected excess returns than the constrainedportfolios in Table 3.4, but also have higher volatility. Note that when γ tends to zero, γ →
0, thisis the log-optimal case lim γ → E h γ ( W γT − (cid:12)(cid:12) W t = w, Z t = z i , but in the formula for the optimal π ∗ we can simply set γ = 0 to obtain the log optimal. Before moving on it is important to pointout that the high Sharpe ratios shown in Tables 3.2, 3.3 and 3.4 are obtained from in-sample back-testing, and for real-life traders there are out-of-sample issues such as parameter error that havesignificant effect on portfolio performance; the data analysis shown here is merely to gain someintuition on statistical arbitrage through the lens of our multiple co-integrated model.Finally, we see in Figure 3.5 the time series of the logarithm of the optimal wealth processes.The case considered is the model with co-integration, for the portfolio with risk aversion coefficient γ = −
1, and the overall proportion of portfolio weight in the risky assets reduced to 0.5% of thatprescribed by the optimal, in other words π ∗ ← . × π ∗ . Note that reducing the risky allocationdoes not affect the Sharpe ratio. From the figure we can observe that the optimal co-integration-based portfolio out-preforms the SPY ETF – even during the 4th quarter of 2015 when there wasthe devaluation of the Renminbi of China. The 2015 devaluation was a series of events executedsurprisingly by the People’s Bank of China, which drove stock markets in the United States, Europe,and Latin America into decline. The time series shown in Figure 3.5 demonstrates that the optimalportfolio under the co-integration model can have some immunity – or even generate profits – whenthe market has this type of downward fluctuation. Statistics / Assets Optimal Wealth γ → γ = − γ = − γ = − γ = − E (cid:16) ∆ W t ∆ tW t − r (cid:17) . . . . . V (cid:16) ∆ W t √ ∆ tW t (cid:17) . . . . . Sharpe Ratio . . . . . Table 3.2:
Merton Portfolio Statistics.
Annualised statistics of the multivariate Merton portfo-lio wealth processes without co-integration, where E (cid:16) ∆ W t ∆ tW t − r (cid:17) is the expectation of excess return, V (cid:16) ∆ W t √ ∆ tW t (cid:17) is the volatility of excess return. tatistics / Assets Optimal Wealth γ → γ = − γ = − γ = − γ = − E (cid:16) ∆ W t ∆ tW t − r (cid:17) . . . . . V (cid:16) ∆ W t √ ∆ tW t (cid:17) . . . . . Sharpe Ratio . . . . . Table 3.3:
Statistical-Arbitrage Portfolio Statistics (no constraints).
Annualised statisticsof the wealth processes for unconstrained model, where E (cid:16) ∆ W t ∆ tW t − r (cid:17) is the expectation of excessreturn, V (cid:16) ∆ W t √ ∆ tW t (cid:17) is the volatility of excess return. Statistics / Assets SPY Optimal Wealth γ = − ETF π ⊤ = 0 π ⊤ s = 0 E (cid:16) ∆ W t ∆ tW t − r (cid:17) . . . V (cid:16) ∆ W t √ ∆ tW t (cid:17) . . . Sharpe Ratio . . . Table 3.4:
Statistical-Arbitrage Portfolio Statistics (with constraints).
Annualised statis-tics of the wealth processes for constrained model, where E (cid:16) ∆ W t ∆ tW t − r (cid:17) is the expectation of excessreturn, V (cid:16) ∆ W t √ ∆ tW t (cid:17) is the volatility of excess return. Log o f O p t i m a l W ea l t h Optimal Wealth versus SPY, = -1
UnconstrainedMarket neutral constraint, T s = 0Dollar neutral constraint, T Figure 3.5: Wealth processes for different optimal portfolios.
We have presented a multi-asset model for co-integration that can be used to analyse statisticalarbitrage opportunities. We compute optimal portfolios for both unconstrained stochastic controlproblem and constrained stochastic control problem. The optimal value functions are the solutionsto Hamilton-Jacobi-Bellman equations, which for power utility function, can be solved with anexponential-affine ansatz that leads to a system of ordinary differential equations. This system24onsists of a matrix Riccati equation and two first-order linear ordinary differential equations. Wepresent the stability analyses of the solutions to these differential equations for the constrained andunconstrained cases. We then apply the optimal formulae to historical data and estimate the modelparameters, and find that stability in fact holds given real-life parameter estimates. We solve theordinary differential equations and look at portfolios based on these historical parameters, fromwhich we draw conclusions about implementation of statistical arbitrage strategies. Our first mainconclusion is that, if a short trading time is allotted, then there is significant probability that thespreads will not converge, and hence the optimal value function will have noticeable marginal gainif more trading time is allotted. Our second conclusion is that there is significant in-sample out-performance of statistical arbitrage portfolios over standard multivariate Merton portfolios, but westress parameter estimation error is likely to change this finding for out-of-sample tests.
References [1] H. Abou-Kandil, G. Freiling, V. Ionescu, and G. Jank.
Matrix Riccati equations in control andsystems theory . Systems & Control: Foundations & Applications. Birkh¨auser Basel, 2003.[2] A. Al-Aradi and S. Jaimungal. Outperformance and tracking: Dynamic asset allocation foractive and passive portfolio management.
Applied Mathematical Finance , 25(3):268–294, 2018.[3] B. Angoshtari.
Stochastic modeling and methods for portfolio management in cointegratedmarkets . PhD thesis, University of Oxford, 2014.[4] M. Avellaneda and J. H. Lee. Statistical arbitrage in the us equities market.
QuantitativeFinance , 10(7):761–782, August 2010.[5] H. Chen, S. Chen, Z. Chen, and F. Li. Empirical investigation of an equity pairs tradingstrategy.
Management Science , Articles in Advance:1–20, September 2017.[6] M. A. Davis and S. Lleo. Risk-sensitive benchmarked asset management.
Quantitative Finance ,8(4):415426, June 2008.[7] M. A. Davis and S. Lleo.
Risk-sensitive investment management , volume 19 of
Advanced Serieson Statistical Science & Applied Probability . World Scientific Publishing, September 2014.[8] R. J. Elliott, J. Van Der Hoek, and W. P. Malcolm. Pairs trading.
Quantitative Finance ,5(5):271–276, June 2005.[9] R. F. Engle and Clive W. J. Granger. Co-integration and error correction representationestimation and testing.
Econometrica , 55(02):251–276, 1987.[10] S. Fallahpour, H. Hakimian, K. Taheri, and E. Ramezanifar. Pairs trading strategy optimiza-tion using the reinforcement learning method: a cointegration approach.
Soft Computing ,20(12):5051–5066, December 2016.[11] W. Fleming and M. Soner.
Controlled Markov processes and viscosity solutions , volume 25 of
Stochastic modelling and applied probability . Springer-Verlag, 2nd edition, 2006.[12] A. Galenko, E. Popova, and I. Popova. Trading in the presence of cointegration.
Journal ofAlternative Investments , 15(1):85–97, 2012.2513] E. G. Gatev, W. N. Goetzmann, and K. G. Rouwenhorst. Pairs trading: performance of arelative average arbitrage rule.
The Review of Financial Studies , 19(3):797–827, October 2006.[14] R. A. Horn and C. R. Johnson.
Topics in Matrix Analysis . Cambridge University Press, NewYork, NY, 1991.[15] S. Lee and A. Papanicolaou. Pairs trading of two assets with uncertainty in co-integration’slevel of mean reversion.
International Journal of Theoretical and Applied Finance , 19(8),November 2016.[16] Y. Lei and J. Xu. Costly arbitrage through pairs trading.
Journal of Economic Dynamics andControl , 56:1–19, July 2015.[17] T. Leung and X. Li. Optimal mean reversion trading with transaction costs and stoploss exit.
International Journal of Theoretical and Applied Finance , 18(3), May 2015.[18] T. N. Li and A. Tourin. Optimal pairs trading with time-varying volatility.
InternationalJournal of Financial Engineering , 03(02), 2016.[19] S. Mudchanatongsuk, J. A. Primbs, and W. Wong. Optimal pairs trading: a stochastic controlapproach. In
Proceedings of 2008 American Control Conference . Institute of Electrical andElectronics Engineers, June 2008.[20] A. Tourin and R. Yan. Dynamic pairs trading using the stochastic control approach.
Journalof Economic Dynamics and Control , 37(10):1947–2156, 2013.[21] G. Vidyamurthy.
Pairs trading quantitative methods and analysis . Wiley Finance. John Wiley& Sons, Inc., Hoboken, New Jersey, U.S.A., August 2004.[22] H. Wielandt. On the eigenvalues of a+b and ab.
Journal of Research of the National Bureauof Standards, Section B: Mathematical Sciences , 77B(1 and 2):61–72, January 1973.[23] W. M. Wonham. On a matrix Riccati equation of stochastic control.