High-dimensional statistical arbitrage with factor models and stochastic control

TL;DR

This paper develops a high-dimensional statistical arbitrage framework combining factor models and stochastic control, providing explicit strategies for market-neutral portfolios with practical constraints and extensive simulations.

Contribution

It introduces a novel analytical approach to high-dimensional arbitrage using factor models and stochastic control, including handling constraints and transaction costs.

Findings

Closed-form optimal strategies for high-dimensional portfolios.

Effective incorporation of constraints like dollar-neutrality.

Successful Monte Carlo simulations with 100 assets.

Abstract

The present paper provides a study of high-dimensional statistical arbitrage that combines factor models with the tools from stochastic control, obtaining closed-form optimal strategies which are both interpretable and computationally implementable in a high-dimensional setting. Our setup is based on a general statistically-constructed factor model with mean-reverting residuals, in which we show how to construct analytically market-neutral portfolios and we analyze the problem of investing optimally in continuous time and finite horizon under exponential and mean-variance utilities. We also extend our model to incorporate constraints on the investor's portfolio like dollar-neutrality and market frictions in the form of temporary quadratic transaction costs, provide extensive Monte Carlo simulations of the previous strategies with 100 assets, and describe further possible extensions of our work.

Figures

Equations

$$d{R}_t=\Lambda d{F}_t+d{X}_t,$$

$$d{F}_t=\tilde{\Lambda }d{R}_t$$

$$d{X}_t=A(\mu -X_t)dt+\sigma d{B}_t,$$

$$d{R}_{ti}=\sum _j{\Lambda }_{ij}d{F}_{tj}+d{X}_{ti},$$

$$d{F}_{tj}=\sum _k{\tilde{\Lambda }}_{jk}d{R}_{tk}.$$

$$d{R}_{ti}=\sum _k\left(\sum _j{\tilde{\Lambda }}_{jk}{\Lambda }_{ij}\right)d{R}_{tk}+d{X}_{ti}=\sum _k{c}_{ik}d{R}_{tk}+d{X}_{ti}.$$

$${\tilde{p}}_i:=(-c_{i1},-c_{i2},\dots ,-c_{i,i-1},1-c_{ii},-c_{i,i+1},\dots ,-c_{iN})$$

$$\sum _i({\pi }_{it}{\tilde{p}}_i)\cdot \text{mathds}1=\sum _i{\pi }_{it}({\tilde{p}}_i\cdot \text{mathds}1)={\pi }_t\cdot p$$

$$d{W}_t={\pi }_t\cdot d{X}_t+(W_t-{\pi }_t\cdot p)rdt$$

$$H(t,x,w)=\underset{\pi \in {\mathcal{A}}_{[t,T]}}{\sup }{\mathbb{E}}_{t,x,w}\left[u(W_T)\right]$$

$$d{W}_t=\left({\pi }_t^{\prime }A(\mu -X_t)+(W_t-{\pi }_t^{\prime }p)r\right)dt+{\pi }_t^{\prime }\sigma d{B}_t$$

$$d{X}_t=A(\mu -X_t)dt+\sigma d{B}_t,$$

$$\begin{gathered} 0={\partial }_{t}H+(\mu -x{)}^{\prime }{A}^{\prime }{\nabla }_{x}H+\frac{1}{2}\mathrm{Tr}(\sigma {\sigma }^{\prime }{\nabla }_{xx}H)+ \\ \underset{\pi }{\sup }\left(\left({\pi }^{\prime }A(\mu -x)+(w-{\pi }^{\prime }p)r\right){\partial }_{w}H+\frac{1}{2}{\pi }^{\prime }\sigma {\sigma }^{\prime }\pi {\partial }_{ww}H+{\pi }^{\prime }\sigma {\sigma }^{\prime }{\nabla }_{xw}H\right) \end{gathered}$$

$${\pi }_t^{\ast }=(\sigma {\sigma }^{\prime }{)}^{-1}\frac{A(\mu -X_t)-pr}{\gamma e^{r(T-t)}}+\frac{A^{\prime }(\sigma {\sigma }^{\prime }{)}^{-1}}{\gamma e^{r(T-t)}}\left((A(\mu -X_t)-pr)(T-t)-Apr\frac{(T-t{)}^2}{2}\right).$$

$$h(t,x)={\mathbb{E}}_{t,x}^{\ast }\left[{\int }_t^T\frac{1}{2\gamma }(A(\mu -Y_s)-pr{)}^{\prime }(\sigma {\sigma }^{\prime }{)}^{-1}(A(\mu -Y_s)-pr)ds\right]$$

$${\pi }^{\ast }=-(\sigma {\sigma }^{\prime }{)}^{-1}\frac{\mathcal{D}H}{{\partial }_{ww}H}$$

$$4\underset{0\le s\le T}{\max }||{\Lambda }_0(s)||<1\text{and}32\underset{0\le s\le T}{\max }||{\Lambda }_1(s)||<1,$$

$$C_0(s)=A^{\prime }(\sigma {\sigma }^{\prime }{)}^{-1}A(I_N+A(T-s)),C_1(s)=A^{\prime }(\sigma {\sigma }^{\prime }{)}^{-1}(I_N+A(T-s))\sigma$$

$$\Omega (s)={\int }_0^s{e}^{-A(s-u)}\sigma {\sigma }^{\prime }{e}^{-A^{\prime }(s-u)}du.$$

$$H(t,x,w)=\underset{\pi \in {\mathcal{A}}_{t,T}}{\sup }{\mathbb{E}}_{t,x,w}\left[W_T-{\int }_t^T\frac{\gamma (s)}{2}\frac{d}{d\tau }{{\mathrm{Var}}_s(W_{\tau })|}_{\tau =s}ds\right]$$

$${\pi }_t^{\ast }=(\gamma (t)\sigma {\sigma }^{\prime }{)}^{-1}\left(A(\mu -X_t)-pr\right)e^{r(T-t)}.$$

$${\pi }_t^{\ast }=(\gamma (t)\sigma {\sigma }^{\prime }+\alpha (t)p{p}^{\prime }{)}^{-1}\left(A(\mu -X_t)-pr\right)e^{r(T-t)}.$$

$${\pi }_t^{\ast }=(\gamma e^{r(T-t)}\sigma {\sigma }^{\prime }+\alpha (t)p{p}^{\prime }{)}^{-1}\left(A(\mu -X_t)-pr-\gamma \sigma {\sigma }^{\prime }(b(t)+c(t)X_t)\right).$$

$$0={\partial }_{t}c-A^{\prime }c-cA-\gamma c\sigma {\sigma }^{\prime }c+e^{r(T-t)}(A+\gamma \sigma {\sigma }^{\prime }c{)}^{\prime }M(t)(A+\gamma \sigma {\sigma }^{\prime }c)$$

$$0={\partial }_{t}b-A^{\prime }b+cA\mu -e^{r(T-t)}(A+\gamma \sigma {\sigma }^{\prime }c)M(t)(A\mu -pr-\gamma \sigma {\sigma }^{\prime }b)-\gamma c\sigma {\sigma }^{\prime }b,$$

$$H(t,x,w,\pi )=\underset{I\in {\mathcal{A}}^{\ast }}{\sup }{\mathbb{E}}_{t,x,w,\pi }\left[W_T-{\int }_t^T\frac{\gamma (s)}{2}\frac{d}{d\tau }{{\mathrm{Var}}_s(W_{\tau })|}_{\tau =s}ds-\frac{1}{2}{\int }_t^T{I}_s^{\prime }C{I}_{s}ds\right].$$

$$I_t=\mathrm{Aim}(t,X_t)+\mathrm{Rate}(t){\pi }_t,$$

$$\mathrm{Aim}(t,x)={\int }_t^{T}f(s){\mathbb{E}}_{t,x}[\mathrm{Frictionless}(s)]ds$$

$${\pi }_s^{\ast }={\pi }_t+{\int }_t^s:\exp \left({\int }_u^s\mathrm{Rate}(v)dv\right):\mathrm{Aim}(u,X_u)du,$$

$$\mathrm{Rate}(t)=C^{-1/2}D\tanh (D(t-T))C^{1/2}$$

Full text

Keywords: statistical arbitrage, factor models, algorithmic trading, Ornstein-Uhlenbeck process, mean reversion, stochastic control.

Word count: 10,318 words.

1 Introduction

Modeling of pairs trading based on stochastic control has been an active research topic in mathematical finance for the last few years. After the papers by Jurek and Yang (2007) and Mudchanatongsuk et al. (2008), an increasing number of models have been proposed in this framework (see, for example, Chiu and Wong (2011), Yan and Tourin (2013), Liu and Timmermann (2013)), in which generally they assume that some statistically-designed relation between the prices of two assets is a mean-reverting stochastic process and find a dynamic optimal allocation in continuous time in some version of the classical Merton framework. More recently, a number of papers have also studied the optimal entry and exit points when trading a couple of cointegrated assets, such as Leung and Li (2015), Lei and Xu (2015), Ngo and Pham (2016), and Kitapbayev and Leung (2017).

In the high-dimensional case, however, relatively little model-based research has been conducted. Cartea and Jaimungal (2016) and Lintilhac and Tourin (2016) investigate a multidimensional generalization of the model in Yan and Tourin (2013) and apply stochastic control to solve a Merton-like problem in continuous time on a collection of cointegrated assets, with exponential utility and finite horizon. In a different direction which is not exactly statistical arbitrage, Cartea et al. (2018) addresses an optimal execution problem on a basket of multiple cointegrated assets, which they also solve with control techniques. Finally, without using stochastic control, the paper by Avellaneda and Lee (2010) carries out a data-based study of statistical arbitrage in the US equity market by proposing a factor model with mean-reverting residuals and a threshold-based strategy. This model is further analyzed and extended by Yeo and Papanicolaou (2017), who discuss risk control and develop an optimization method to allocate the investments given the trading signals.

The previous papers in this high-dimensional framework thus either apply stochastic control to a given mean-reverting process or use a factor model to construct this process and then choose the trading signals based on threshold rules, but none of them considers the combination of these two techniques. The present paper aims to fill this gap by providing a study of statistical arbitrage in a high-dimensional setting that combines factor models and the tools from stochastic control, considers transaction costs and statistical arbitrage constraints, and obtains closed-form optimal strategies which are interpretable and easy to implement computationally.

More precisely, in our framework an investor observes the returns of a high-dimensional collection of risky assets and, similar to Avellaneda and Lee (2010) and Yeo and Papanicolaou (2017), uses historical data to statistically construct a factor model such that the cumulative residuals are assumed to be mean-reverting and following an Ornstein-Uhlenbeck process. However, unlike these previous studies, these residuals may be correlated and interdependent and, based on their behavior, the investor must decide how to optimally allocate her wealth in the risky assets and a riskless security so that the expected utility of her terminal wealth is maximized and she is market-neutral222In this paper we use the expression “market-neutral” as in Avellaneda and Lee (2010) to refer to factor neutrality.. There are three main results in this paper:

First, for a big class of statistically-constructed factor models that includes PCA we show how the investor may theoretically construct market-neutral portfolios without solving any optimization problem (unlike the approach followed in Yeo and Papanicolaou (2017) or Boyd et al. (2017), for example) provided that the factor model holds, and we show how this makes the optimal allocation problem analytically tractable and guarantees market-neutrality by construction. These portfolios are explicitly computable and depend quadratically on the factor model loadings and, to the best of our knowledge, using this construction to connect factor models and stochastic control theory is new.

Second, using these explicit market-neutral portfolios as control variables, we show how the investor should trade optimally in continuous time to maximize either an exponential utility or a mean-variance objective, obtaining explicit analytic forms of the optimal strategies in both cases in this high-dimensional setting. The structure of these optimal strategies is related to the classical solution of the Merton problem and is affine in the deviation of the residuals from their statistical mean, thus giving a precise estimate of how much we should buy when the assets are underpriced and how much we should sell when they are overpriced, as in classical pairs trading. The coefficients are given by the solution of matrix Riccati differential equations and depend quadratically on the factor model loadings, and the strategies in both the exponential and the mean-variance case are surprisingly similar except for a non-myopic correction term that does not appear in the classical framework under a geometric Brownian motion. This arises from the fact that in our case the drift of the underlying Ornstein-Uhlenbeck process is stochastic.

The structure and the techniques to find these affine strategies are thus similar in spirit to those in the affine process literature in finance (see Duffie et al. (2003) for a broad survey), to the more recent affine control literature in algorithmic trading (see, for example, Cartea et al. (2015) and the references therein), and to the literature on extensions of the Merton problem (see, for example, Benth and Karlsen (2005), Liang et al. (2011), Fouque et al. (2015) and Moutari et al. (2017), which deal with a single risky asset in the context of the Schwarz model or in geometric Brownian motion with stochastic drift or volatility; and Brendle (2006) and Bismuth et al. (2019), which consider the multiasset case in the setting of geometric Brownian motion with uncertain drift). While the techniques that we use to find the optimal strategies are therefore classical, the framework and the results are new because the mean-reverting behavior of the underlying stochastic process arises from the residuals of a factor model and in the context of statistical arbitrage, and we consider the general case of an arbitrary number of assets with a market-neutrality restriction and a general matrix Ornstein-Uhlenbeck process. Moreover, the explicit solutions allow us to understand the dependence of the optimal strategies on specific elements of a statistical arbitrage strategy (such as the factor model, its loadings matrix and its connection with market-neutrality, and the mean-reversion speed of the residuals and their correlation structure), and to compare arbitrageurs with exponential and mean-variance utilities.

Third, we extend the previous results in two directions by discussing how to incorporate into the model soft constraints frequently imposed by arbitrageurs such as dollar-neutrality, and also market frictions in the form of quadratic transaction costs, inspired by Garleanu and Pedersen (2013) and Garleanu and Pedersen (2016) and also by the more general quadratic transaction cost and linear price impact literature in portfolio theory (see, for example, Moreau et al. (2017) and Muhle-Karbe et al. (2017) for some new research directions, and Obizhaeva and Wang (2013), Rogers and Singh (2010), Almgren and Chriss (2001) and Bertsimas and Lo (1998) for some classical papers). In both extensions, we again find explicit analytic strategies which are easily interpretable, and which quantitatively correspond to quadratic corrections in the structure of the original optimal strategies (when adding soft constraints like dollar neutrality) or to “tracking” averages of the future original optimal portfolios (when adding quadratic transactions costs). Moreover, in both cases these new strategies depend quadratically on the loadings of the factor model. Again, the novelty of the results comes from the study of these questions (dollar neutrality, transaction costs, etc.) in a new context in which they are crucial (statistical arbitrage with an arbitrary number of assets following a general matrix Ornstein-Uhlenbeck process and a market-neutrality restriction, in particular using control techniques and a factor model), and this framework and the strategies that we find are new to the best of our knowledge.

To conclude the paper with a more empirical analysis, we also perform extensive numerical simulations with a high-dimensional number of assets. This gives further insights about the behavior of the previous strategies that are not obvious when looking at the corresponding equations, and allows us to understand the sensitivity of the model parameters and the dependence on the underlying factor model. This high-dimensional numerical study is also new with respect to the existing literature, and the main conclusions are that (1) the exponential-utility strategies are more profitable than the mean-variance strategies and they also take more extreme positions, (2) after some initial up and downs the sample paths of the different wealth processes progressively stabilize due to the asymptotic properties of the Ornstein-Uhlenbeck process, (3) increasing the risk-control parameters consistently produces a concentration of the distribution of the terminal wealth around smaller values, and (4) imposing market neutrality when the loadings of the factor model get bigger leads to more aggressive strategies whose terminal wealth has a higher variance.

The remainder of the article is organized as follows. In section 2 we introduce our model, construct the market-neutral portfolios that make the problem analytically tractable, and formulate the control problems. Next, in section 3 we present the basic results under the exponential and the mean-variance frameworks, whereas section 4 extends these results by considering the addition of soft constraints and of quadratic transaction costs. Section 5 contains Monte Carlo simulations that provide further insight about the qualitative behavior of the found strategies, and section 6 presents the main conclusions and proposes future new directions of research. Finally, an appendix contains all the proofs.

2 The model

2.1 Set-up and assumptions

In the remainder of this paper we will consider the following general framework. We will assume that an investor observes the returns of a large number $N$ of risky assets and, like in classical portfolio theory based on stochastic control, she must decide how to dynamically allocate her wealth by investing in them or in a riskless asset with constant interest rate $r$ so that the expected utility of her wealth at a finite terminal time $T$ is maximized. However, unlike the classical framework and the existing literature, to do so she will execute a statistical arbitrage strategy based on a factor model, in which instead of trading depending on the state of the original returns she will trade depending on the behavior of the residuals, which will be the trading signals. For example, in the case of two assets, this is equivalent to classical pairs trading, in which the investor may perform a simple linear regression on the returns of two historically correlated securities and, depending on how far the oscillation of the residual is from its historical average, she decides if there is a mispricing and opens and closes long and short positions in the original assets in a market-neutral way. In this paper, we will study the generalization of this to the high-dimensional case of an arbitrary number of assets, in which we substitute the simple linear regression by a statistical factor model and we study the optimal allocations under the framework of stochastic control, assuming a mean-reverting stochastic model for the behavior of the residuals.

More precisely, we make the following three general assumptions on how the investor will generate these residuals and what dynamics they will have:

1. (1)

Assumption 1: The investor has computed a factor model for the returns of the risky assets, which will hold during the investment finite horizon and is given in differential form by

[TABLE]

where $R_{t}$ is the cumulative asset return process, $\Lambda$ is the (constant-in-time333We make this assumption for analytic tractability in our control framework, and given that the trading frequency is in general higher than the frequency at which these loadings will change significantly.) loadings matrix, $F_{t}$ is the cumulative factor returns process, and $X_{t}$ is the cumulative residual returns process444Here we have written the factor model in a somewhat unusual differential form in terms of the cumulative residuals and returns because of notational simplicity for this section of the paper. In practice, however, the factor model will be estimated in discrete time, by replacing the differentials by the corresponding discrete increments (so, for instance, $dR_{t}$ should be replaced by the daily, hourly, etc. asset returns, $dF_{t}$ would be just the corresponding daily, hourly, etc. factors returns, and so forth). In any case we will only use this notation and framework in this section of the paper, and the reader may look at Avellaneda and Lee (2010) for essentially the same continuous/discrete time framework and some estimation techniques. . 2. (2)

Assumption 2: This factor model has been computed statistically by using some version of PCA555See Lettau and Pelger (2020) for some new versions of PCA that might be particularly interesting for this problem., so the rows of $\Lambda$ are the largest eigenvectors of some square matrix and the discrete-time version of $dF_{t}$ (i.e., the daily, hourly, etc. factors returns) is computed by linearly regressing the discrete-time version of $dR_{t}$ (i.e., the daily, hourly, etc. assets returns) on some rescaling of $\Lambda$, so

[TABLE]

for some matrix $\tilde{\Lambda}$. However, the only fact we will need about this assumption is that (2.2) holds for some matrix $\tilde{\Lambda}$, which allows for a bigger class of factor models than classical PCA. 3. (3)

Assumption 3: The process $X_{t}$ given by the cumulative residuals is mean-reverting. In particular, for analytic tractability we assume that it is a matrix $N$-dimensional Ornstein-Uhlenbeck process satisfying the following stochastic differential equation with known parameters

[TABLE]

where $A$ is a constant $N$-dimensional square matrix, $\mu$ is a constant $N$-dimensional vector, $B_{t}$ is a vector of $m$ independent Brownian motions in the usual complete filtered probability space $(\Omega,\mathcal{F},{\mathbb{P}},(\mathcal{F}_{t})_{0\leq t\leq T})$, and $\sigma$ is a constant $N\times m$ matrix such that the instantaneous covariance matrix $\sigma\sigma^{\prime}$ is invertible.

The previous framework thus combines high-dimensional statistical arbitrage, factor models and stochastic control in a way which is new to the best of our knowledge, and it extends several models in the existing literature. For example, statistical arbitrage models based on a more particular case of Assumptions 1, 2, 3 (in which the residuals are assumed to be independent one-dimensional Ornstein-Uhlenbeck processes, so $A$ and $\sigma$ are diagonal) and in which no stochastic control methods are applied have been studied empirically in the US equity market by Avellaneda and Lee (2010) and Yeo and Papanicolaou (2017). In a different direction, if we consider the particular case of removing the factor model by making $\Lambda=0$, we have the situation in which the returns are globally mean-reverting following a matrix Ornstein-Uhlenbeck process, which has also been studied empirically and analytically using stochastic control techniques in the context of optimal execution in Cartea et al. (2018).

2.2 Making the model tractable and imposing market-neutrality

Unlike the classical literature on portfolio choice based on stochastic control, choosing as control variables the amount of capital that the agent invests in each of the $N$ risky assets of the previous framework might make the optimal allocation problem analytically intractable. Indeed, since we only have information about the dynamics of the residuals and not directly about the returns like in the classical framework, these residuals are not independent, and the factors depend on the returns, the classical approach would lead to complicated interdependencies. Moreover, since the investor is executing a statistical arbitrage strategy, we would need to incorporate additional market neutrality constraints666Here we use market-neutrality in the sense of Avellaneda and Lee (2010). so that the returns of the strategy do not depend on the model factors, but just on the idiosyncratic component of the model given by the residuals. This would complicate the problem further, and would require numerical optimization methods as done in Yeo and Papanicolaou (2017).

In this paper, on the contrary, we approach both problems simultaneously and we solve them analytically by following a new approach. This is based on the following proposition, which shows that, by using the original $N$ risky assets, it is possible to construct analytically $N$ market-neutral portfolios whose returns only depend on one coordinate of $X$, which simplifies the complexity of the problem and makes it analytically tractable:

Proposition 2.1**.**

Under the previous assumptions, it is possible to construct explicitly $N$ market-neutral portfolios such that investing any real number $\pi_{it}$ of dollars in the $i$-th one at time $t$ yields an instantaneous return of $\pi_{it}dX_{ti}$ (and hence a combined return of $\pi_{t}\cdot dX_{t}$).

Moreover, the total amount of capital invested at time $t$ by doing so is $\pi_{t}\cdot p$ for an explicit constant-in-time vector $p\in{\mathbb{R}}^{N}$, which depends quadratically on the factor model loadings.

Proof.

The mathematical construction of the market-neutral portfolios under the given assumptions is surprisingly straightforward and involves just a linear projection. Indeed, (2.1) implies that

[TABLE]

whereas (2.2) yields

[TABLE]

Combining the two previous equations we find that, for $c_{ik}:=\sum_{j}\tilde{\Lambda}_{jk}\Lambda_{ij}$,

[TABLE]

Thus, if at time $t$ we hold the (explicitly constructible) constant-in-time portfolio given by

[TABLE]

(i.e., we invest $-c_{i1}$ dollars in the first asset, $-c_{i2}$ dollars in the second one, and so on), we automatically obtain an instantaneous return of $dX_{ti}$, which is market neutral and depends only on the $i$th coordinate of the process $X_{t}$. Further, from the above equations it is also obvious that for any real number $\pi_{it}$, $\pi_{it}\tilde{p}_{i}$ will also be market-neutral and yielding a return of $\pi_{it}dX_{ti}$, and the same applies to $\sum_{i}\pi_{it}\tilde{p}_{i}$, which will have a return of $\sum_{i}\pi_{it}dX_{ti}=\pi_{t}\cdot dX_{t}$.

Finally, regarding the last part of the statement just observe that the total amount of capital invested in the strategy $\pi_{t}=(\pi_{it})_{1\leq i\leq N}$ at time $t$ is simply

[TABLE]

where $p:=(\tilde{p}_{i}\cdot\mathds{1})_{1\leq i\leq N}$, which concludes our proof. ∎

Remark 2.1**.**

Note in particular that, if $\Lambda$ or $\tilde{\Lambda}$ are sparse matrices, then most of the $c_{ik}$ in the above construction will be 0, so the investor will be investing in a few number of assets in each market-neutral portfolio and this could significantly reduce his transaction costs while rebalancing his positions. In particular, Pelger and Xiong (2018) discusses a way of obtaining this kind of sparse factor model.

The key consequence of the above proposition is that, if we choose as control variables the amount of capital $\pi_{t}$ that we wish to invest in these $N$ market-neutral portfolios (instead of directly in the original assets) at time $t$, the dynamics of the problem get simpler, they only depend separately on the coordinates of $X$, and we have market-neutrality by construction. This solves simultaneously the two problems we discussed before and allows us to connect stochastic control and the factor model in a simple way, and it is therefore the approach which we will adopt in the remainder of this paper.

Note also that, under these new control variables, all the information about the factor model and in particular about its loadings matrix is now encoded in the vector $p$, which will play an important role in the remaining sections. Moreover, some statements about the strategies must be rewritten in terms of it within this new framework. For instance, in the new setting a strategy $(\pi_{t})_{0\leq t\leq T}$ is dollar-neutral at $t$ if $p\cdot\pi_{t}=0$, since as we mentioned before $p\cdot\pi_{t}$ is the total capital spent at time $t$.

2.3 Formulation of the control problems

Under the previous framework, now we formulate rigorously the control problems we will study in the paper. We suppose that the investor executes the following trading strategy: at each time $t\in[0,T]$, she invests $\pi_{t}$ dollars in the risky market-neutral portfolios we constructed in Proposition 2.1, and she invests her remaining wealth in the risk-free asset with constant interest rate $r$, so that the resulting strategy is self-financing. Thus, assuming for the moment no market frictions or constraints (which will be both considered in section 4), the evolution of her wealth is given by the equation

[TABLE]

and she aims to choose $\pi_{t}$ to maximize the expected utility of her terminal wealth, which is given by $u(W_{T})$ for a fixed utility function $u$.

Supposing further that she trades continuously in time, this means that mathematically she must solve the high-dimensional non-linear stochastic optimization problem given by

[TABLE]

subject to

[TABLE]

[TABLE]

where ′ indicates transposition, and the admissible set $\mathcal{A}_{[t,T]}$ is the set of all the $\mathcal{F}_{s}$-predictable and adapted processes $(\pi_{s})_{s\in[t,T]}$ in ${\mathbb{R}}^{N}$ with the minimal technical restrictions that ${\mathbb{E}}[\int_{t}^{T}||\pi_{s}||^{2}ds]<\infty$ (so Ito’s formula may be applied and doubling strategies are excluded) and the above SDEs have a unique strong solution.

Finally, the associated dynamic programming equation of the problem is non-linear and $(N+2)$-dimensional, and is given by

[TABLE]

with terminal condition $H(T,x,w)=u(w)$.

The problem is therefore formally related to the classical Merton framework, but instead of a geometric Brownian motion there is a multidimensional Ornstein-Uhlenbeck process which makes it impossible to combine the dynamics of $W$ and $X$ into a single equation and to get rid of the $N$-dimensional state variable $x$. Moreover, unlike the previous studies on extensions of the Merton problem with an Ornstein-Uhlenbeck process discussed in section 1, in (2.4) and (2.5) the mean-reverting behavior of the underlying stochastic process arises in the context of statistical arbitrage and from the residuals of a factor model (which is encoded in the vector $p$ of the equations above and which will play an important role in the following sections), and we consider the general case of an arbitrary number of assets with a market neutral restriction. Furthermore, the model will be extended in section 4 to incorporate other important features of statistical arbitrage strategies, like dollar neutrality restrictions and transaction costs, and we will analyze the impact of the factor model on these extensions.

Quite surprisingly, the previous problems admit interpretable closed-form solutions – which is computationally useful in this high-dimensional setting, and which allows us to understand the influence of the model parameters and especially of the factor model – in the cases in which the utility is exponential or of a Markowitz-inspired mean-variance type, but not for other usual choices of utility functions, like the HARA family. We show this in the following two sections, first for the simple setup of (2.4) and (2.5) in section 3, and then extending the model in section 4 to incorporate soft constraints on the investor’s portfolio and quadratic transaction costs.

3 The frictionless results

In this section we present the closed-form, optimal strategies for the problem given by (2.4) and (2.5) in the cases in which the utility is exponential or of a mean-variance type, discussing the former in the first subsection and the latter in the second one.

3.1 The exponential utility case

In the first setting, the explicit description of the optimal strategy is given by the following main theorem (see Cartea and Jaimungal (2016) and Lintilhac and Tourin (2016) for related results with an exponential utility):

Theorem 3.1**.**

Under an exponential utility (so $u(w)=-e^{-\gamma w}$ for some $\gamma>0$) and the technical condition described in our verification theorem (Proposition 3.2 below), the optimal portfolio to have at time $t$ according to (2.4) is explicitly computable and given by

[TABLE]

The result follows from the following two propositions, whose proof is given in appendix A.1 using classical stochastic control techniques:

Proposition 3.1** (Solving the PDE).**

The value function $H$ associated to (2.4) and (2.5) when $u(w)=-e^{-\gamma w}$ is explicitly computable and admits the probabilistic representation $H(t,x,w)=-\exp(-\gamma(we^{r(T-t)}+h(t,x)))$ where

[TABLE]

and $dY_{t}=rpdt+\sigma dB_{t}^{*}$ for a new Brownian motion $B^{*}$ under a new probability law ${\mathbb{P}}^{*}$. The associated optimal control in feedback form is then

[TABLE]

*where $\mathcal{D}H=\left(A(\mu-x)-pr\right)\partial_{w}H+\sigma\sigma^{\prime}\nabla_{xw}H$.

Proposition 3.2** (Verification).**

The strategy given in Theorem 3.1. is indeed optimal if

[TABLE]

where $\Lambda_{0}(s)$ and $\Lambda_{1}(s)$ are the diagonal matrices containing, respectively, the eigenvalues of $\Omega^{1/2}(C_{0}+C_{0}^{\prime})\Omega^{1/2}(s)$ and $\Omega^{1/2}C_{1}C_{1}^{\prime}\Omega^{1/2}(s)$, for

[TABLE]

[TABLE]

Besides being a closed-form strategy which is easily implementable in our high-dimensional setting, the above optimal portfolio is also interpretable. Indeed, the first term of the optimal policy is Merton-like since it represents the drift of the underlying process (which here is stochastic unlike in the classical geometric Brownian motion) minus the adjusted risk-free rate (which here depends on the loadings of the factor model via $p$). This is divided by a measure of the volatility (which is given by $\sigma\sigma^{\prime}$, the instantaneous quadratic covariation of $X$) and the Arrow-Pratt coefficient of absolute risk-aversion of the value function with respect to the wealth $w$ (i.e., $-\partial_{ww}H/\partial_{w}H$), which is the product $\gamma e^{r(T-t)}$, where $\gamma$ is the risk aversion parameter of the utility function and the factor $e^{r(T-t)}$ measures the gains from interest between $t$ and $T$.

On the other hand, the second summand is a non-myopic correction term which again depends linearly on the drift of $X$, and whose effect vanishes when we approach the terminal time $T$. Moreover, while the first term does not depend explicitly on the terminal time $T$, this correction term does, reflecting the fact that, since there are non-zero interest rates and moreover the behavior of the residuals is oscillating, the investor must keep in mind the final horizon to decide if she bets on the mean-reversion cycle before that time. Finally, as the risk-aversion parameter $\gamma$, the instantaneous volatility $\sigma\sigma^{\prime}$, or the interest rate $r$ increase, the optimal portfolio vector $\pi_{t}^{*}$ gets closer to 0, implying that the investor will simply invest most of her wealth in the riskless asset.

The above strategy is also intuitive within our framework of statistical arbitrage with a factor model, and sheds further light on the problem. Indeed, note that the current state of the residual process $X_{t}$ only appears in the strategy through the terms in $A(\mu-X_{t})$, which essentially tells us to invest more in the risky assets the further their residuals are from their historical mean $\mu$ and in a way proportional to the historical mean reversion speed given by $A$, like in classical pairs trading. Moreover, all the remaining terms depend jointly on the factor model and the interest rate through the term $pr$, which reflects the cost of the leverage associated to imposing market-neutrality through the loadings of the factor model. In particular, note that, the bigger the loadings of the factor model are and hence the bigger $p$ is, the more the agent will need to invest to achieve market neutrality (again like in pairs trading with a big beta) and the bigger her leverage will be, and this will affect the optimal strategy depending on the interest rate $r$.

Finally, regarding the technical optimality conditions, intuitively they arise from the fact that $H(t,X_{t},W_{t}^{*})$, the value function evaluated at the wealth process $W_{t}^{*}$ corresponding to the optimal strategy, may blow up because of the exponential function coming from the exponential utility. In particular, since $W_{t}^{*}$ ends up being an Ito process depending quadratically on $X_{t}$ and $X_{t}$ is Gaussian, the term $\exp(-\gamma W_{t}^{*}e^{r(T-t)})$ is related to the moment generating function of a chi-squared distribution, which blows up far away from 0. Thus, these technical conditions are just ensuring that the corresponding functions are integrable. Interestingly, this does not depend on the risk-aversion parameter $\gamma$, the interest rate $r$, or the factor model used (captured by $p$), but just on the parameters of $X$ and the terminal time $T$.

3.2 The mean-variance case

In the second, Markowitz-inspired mean-variance framework, the investor tries to maximize her expected terminal wealth, but at the same time she continuously penalizes at each instant the instantaneous variance of her wealth process according to a time-dependent volatility-aversion function $\gamma(t)$. The optimal strategy in this case is again available in closed form and interpretable and, for an appropriate choice of this volatility-aversion function, we obtain the same optimal policy as in the exponential case but without the correction term. This is shown in the following theorem, whose proof is given in appendix A.2:

Theorem 3.2**.**

If $\gamma(t)$ is continuous and positive on $[0,T]$, the problem in (2.4) with the following mean-variance objective function

[TABLE]

has explicit optimal portfolio at $t$ given by

[TABLE]

In particular, for $\gamma(t)=\gamma e^{2r(T-t)}$, the above optimal policy is the same as the first term of the corresponding one in Theorem 3.1.

Regarding the interpretation of the mean-variance strategy within our context of statistical arbitrage and its connection with the exponential-utility arbitrageur, there are two important remarks.

First, as we mentioned, the new optimal strategy is the same as the myopic part of the exponential case modulo the value of $\gamma(t)$. In particular, this means that, unlike the exponential arbitrageur, the mean-variance arbitrageur will not take into account the expected number of mean-reversion cycles until the terminal time $T$. Moreover, for a non-zero interest rate and a constant volatility aversion $\gamma(t)$, the mean-variance arbitrageur is more aggressive than the corresponding exponential investor with the same $\gamma$, since she will invest significantly more capital (quantitatively, by a factor of $e^{2r(T-t)}$) in going long or short, taking more aggresive positions the higher the interest rate is and the sooner it is with respect to the terminal date.

Second, the optimal strategy has two components like in section 3.1: one term in $A(\mu-X_{t})$ which measures how far the residuals are from their historical mean and how fast they will mean-revert (like in classical pairs trading), and a second term in $pr$ linked to the factor model, which measures the cost of the leverage associated to imposing market neutrality. In particular, note that, the bigger the loadings of the factor model are (and hence the bigger $p$ is), the more aggressive the positions will be and the more leverage the investor will have if $r\neq 0$.

4 Two extensions

In this section of the paper, we complete the picture described in the previous two sections by considering two important extensions within the context of statistical arbitrage with a factor model. In the first subsection, we show how to incorporate in the above strategies soft constraints frequently imposed by arbitrageurs with the example of dollar-neutrality. In the second one, we introduce market frictions in the form of quadratic transaction costs. In both cases, we obtain again closed-form analytic solutions which are interpretable, convenient from a computational perspective in our high-dimensional setting, and which shed further light on the influence of the factor model and its connection with market neutrality.

4.1 Incorporating soft constraints on the admissible portfolios

While imposing restrictions on the portfolios by introducing hard constraints directly on the admissible set $\mathcal{A}_{t,T}$ leads in general to control problems that must be solved numerically (and hence potentially unfeasible in a high-dimensional setting since in applications the number of assets ranges from hundreds to thousands), it is still possible to impose many additional soft constraints in the two frameworks of section 3 without significantly increasing the difficulty of the problems, by just adding a carefully chosen penalty term to the corresponding objective function.

As an illustration of this, we give in the next corollary the corresponding optimal strategies when a dollar-neutrality restriction is softly enforced. To do so, recall that, within the framework of section 2 that imposed market-neutrality within the factor model, a strategy $\pi_{t}$ is dollar neutral if $p\cdot\pi_{t}=0$, which means that the total amount of capital invested at time $t$ is 0. Hence, we can softly enforce dollar neutrality by replacing the wealth process of Theorems 3.1 and 3.2 by the penalized wealth process defined by $d\tilde{W}_{t}:=dW_{t}-\alpha(t)(\pi_{t}\cdot p)^{2}/2dt$ for a certain general time-dependent penalty function $\alpha(t)$. This penalizes non dollar-neutrality (i.e., $\pi_{t}\cdot p\neq 0$) at each time and is quadratic to simplify the optimization process.

The proof follows the same lines as in the previous two cases and is obtained from them by small modifications, so we will omit it for the sake of brevity.

Corollary 4.1**.**

Suppose that dollar neutrality is softly enforced by replacing the wealth process of Theorems 2.3.1 and 2.3.2 by the penalized wealth process defined by $d\tilde{W}_{t}:=dW_{t}-\alpha(t)(\pi_{t}\cdot p)^{2}/2dt$. Then

1. (1)

The problem with mean-variance utility has optimal portfolio at $t$ given by

[TABLE] 2. (2)

The problem with exponential utility has optimal portfolio at $t$ given by

[TABLE]

where $c(t)$ is an $N\times N$ symmetric matrix and $b(t)$ is an $N$-dimensional vector, vanishing when $t\rightarrow T$, and with coordinates depending on $A,\sigma,rp,\gamma$ and $\alpha(t)$. In particular, $c(t)$ is given by the solution of the matrix Riccati ODE

[TABLE]

and $b(t)$ is the solution of the linear system of ODEs

[TABLE]

both with terminal conditions $b(T)=c(T)=0$ and where $M(t)=(\gamma\sigma\sigma^{\prime}e^{r(T-t)}+\alpha(t)pp^{\prime})^{-1}$.

The resulting optimal policies have therefore the same structure as the two previous strategies of section 3, but now the additional term $\alpha(t)pp^{\prime}$ has been introduced in the inverse to enforce the dollar-neutrality condition. This again depends on the factor model via $p$ and is related to how extreme the capital positions will be because of the market-neutrality restriction, which depends directly on the loadings matrix and hence on $p$. Note in particular that, the bigger the loadings of the factor model are, the bigger $\alpha(t)pp^{\prime}$ will be and hence the bigger the impact of the dollar neutrality restriction will be.

4.2 Incorporating quadratic transaction costs

In this subsection, we extend our model to incorporate market frictions in the form of transaction costs, which play a crucial role when executing statistical arbitrage strategies. We consider in particular quadratic transaction costs, which are in general a measure of price impact or illiquidity and which make the model anaytically tractable.

To do so, rather than looking at the amount of capital $\pi_{t}$ invested in the risky assets at time $t$ as the control variables, we consider the trading intensity $I_{t}$ at which these investements will be made at time $t$, which is given by $d\pi_{t}=I_{t}dt$. We can now adapt to our setting the model for transaction costs introduced in Garleanu and Pedersen (2016), who posit, providing market microstructural justification and referring to empirical research, that these transaction costs at time $t$ may be represented quadratically as $I_{t}^{\prime}CI_{t}$ for a certain symmetric positive-definite matrix $C$777The assumption that $C$ is symmetric is without loss of generality, since if the transaction costs are given by $I_{t}^{\prime}\tilde{C}I_{t}$ for a non-symmetric $\tilde{C}$, then one can see that by considering the symmetric part of $\tilde{C}$ (given by $C:=(\tilde{C}+\tilde{C}^{\prime})/2$) we have that $I_{t}^{\prime}\tilde{C}I_{t}=I_{t}^{\prime}CI_{t}$., which essentially comes from the assumption that the price impact of the investor’s actions is linear on its trading intensity $I_{t}$.

Under this framework, we rewrite the performance criteria of Theorem 3.2 by incorporating the adverse effect caused by these transaction costs on the investor’s wealth as a running penalty, obtaining the stochastic optimization problem given by

[TABLE]

As we mentioned, in this new problem the control variable is $I$; $t,x,w,\pi$ are now state variables; and $\mathcal{A}^{*}$ is the set of all $\mathcal{F}$-adapted predictable processes $I_{t}$ such that the corresponding SDEs have a unique strong solution for any initial data and both $I_{t}$ and the resulting $\pi_{t}$ given by $d\pi_{t}=I_{t}dt$ are again in $L^{2}(\Omega\times[0,T])$. Thus, the investor aims to maximize her terminal wealth, but penalizing at each instant both for the risk of her strategy (measured by the instantaneous variance of her wealth process) and for the price impact caused by her actions (reflected in the quadratic transaction costs).

In this new setting, it is again possible to find explicitly the optimal strategy that the investor should follow, which is described in detail in the next theorem:

Theorem 4.1**.**

If $\gamma(t)\geq 0$ and is continuous, the optimal strategy in the above problem “tracks” a moving aim portfolio $\mathrm{Aim}(t,X_{t})$ with a tracking speed of $\mathrm{Rate}(t)$, according to the following ODE describing the evolution of the optimal trading intensity $I_{t}=d\pi_{t}/dt$

[TABLE]

where $\mathrm{Rate}(t)$ is a $N\times N$ negative-definite matrix tending to 0 when $t\rightarrow T$888and given by the solution of a matrix Riccati ODE specified in the Porposition 2.4.2 below., and $\mathrm{Aim}(t,X_{t})$ admits the probabilistic representation

[TABLE]

where $\mathrm{Frictionless}(s)$ is the optimal portfolio at time $s$ in the frictionless case of section 3.2. and $f(s)$ is a certain averaging function given in Proposition 4.3 below.

Furthermore, the optimal portfolio is then

[TABLE]

*where the notation $:\exp(\int_{u}^{t}\cdot ds):$ represents the time-ordered exponential999Recall that the time-ordered exponential of a time-dependent matrix $A(s)$ is defined as $:\exp(\int_{u}^{t}A(s)ds):=\lim_{||\mathcal{P}||\downarrow 0}\prod_{i=1}^{n_{\mathcal{P}}}\exp(A(t_{i})\Delta t_{i})$, where $\mathcal{P}:=\{u=t_{0},t_{1},\ldots,t_{n}=t\}$ is a partition of $[u,t]$, $\Delta t_{i}:=t_{i}-t_{i-1}$, and the product is ordered increasingly in time. If $A(s)$ is a scalar, then obviously $:\exp(\int_{u}^{t}A(s)ds):=\exp(\int_{u}^{t}A(s)ds)$..

Remark 4.1**.**

If in particular the investor has constant volatility aversion (so $\gamma(t)=\gamma$), the matrix Riccati ODE is explicitly solvable and

[TABLE]

for $D:=(\gamma C^{-1/2}\sigma\sigma^{\prime}C^{-1/2})^{1/2}$. Moreover, if the transaction costs are proportional to the volatility (i.e., $C=\lambda\sigma\sigma^{\prime}$ for $\lambda>0$, see Garleanu and Pedersen (2013, 2016) for a market microstructural justification) then this rate is indeed a scalar given by $\sqrt{\frac{\gamma}{\lambda}}\tanh\left(\sqrt{\frac{\gamma}{\lambda}}(t-T)\right)$ and $:\exp\left(\int_{u}^{s}\mathrm{Rate}(v)dv\right):=\cosh\left(\sqrt{\frac{\gamma}{\lambda}}(s-T)\right)/\cosh\left(\sqrt{\frac{\gamma}{\lambda}}(u-T)\right)$.

The result follows from the following sequence of three propositions, which are proved in appendix A.3:

Proposition 4.1** (Conjectured solution).**

The solution of the HJB equation associated to the problem is $H(t,x,w,\pi)=e^{r(T-t)}w+\frac{1}{2}\pi^{\prime}a(t)\pi+\pi^{\prime}b(t,x)+d(t,x)$ if there exist a $N\times N$ symmetric matrix $a(t)$, a $N$-dimensional vector $b(t,x)$ and a scalar function $d(t,x)$ satisfying

1. (1)

The matrix Riccati ODE

[TABLE]

with terminal condition $a(T)=0$. 2. (2)

The vector-valued and the scalar linear parabolic PDEs

[TABLE]

[TABLE]

with terminal conditions $b(T,x)=d(T,x)=0$ and where $\mathcal{L}_{X}$ is the infinitesimal generator of $X$, acting coordinatewise.

*The hypothesized optimal trading intensity at $(t,x,w,\pi)$ is then $I^{*}=C^{-1}(a(t)\pi+b(t,x))$.

Proposition 4.2** (Existence of solutions).**

.

1. (1)

If the volatility aversion $\gamma(t)\geq 0$ and is continuous, then the Riccati equation (4.2) has a symmetric, bounded and negative definite solution on all $[0,T]$. In particular, for $\gamma(t)=\gamma$, the solution is

[TABLE]

for $D:=(\gamma C^{-1/2}\sigma\sigma^{\prime}C^{-1/2})^{1/2}$. 2. (2)

Moreover, under this condition the parabolic PDEs (4.3) and (4.4) have a unique solution satisfying a polynomial growth condition in $x$, and this solution admits the probabilistic representation

[TABLE]

[TABLE]

Furthermore, $b$ has linear growth in $x$ and $d$ has quadratic growth in $x$, both uniformly in $t$.

Proposition 4.3** (Verification).**

Under the assumptions of the previous proposition, the trading intensity given in Theorem 4.1 is indeed optimal with the choices

[TABLE]

The interpretation of the above strategy, which is again explicit and hence easily implementable in practice, is intuitive and complementary to the infinite-horizon model of Garleanu and Pedersen (2016): in the presence of quadratic transactions costs, the investor trades with a certain decreasing rate $\mathrm{Rate}(t)$ towards an aim portfolio $\mathrm{Aim}(t,X_{t})$ depending on the time and the mean-reversion state of the signals $X_{t}$. This aim portfolio is given by a weighted average of the future optimal strategies in the frictionless case, reflecting the fact that now trading is not free and thus to enter a trade the investor must weight the future outcomes derived from the strategy. Moreover, as shown in the above remark, the trading rate is bounded by 1 because of the properties of $\tanh$, depends on $t$ unlike the infinite-horizon case, and is naturally decreasing in the transaction cost parameter $\lambda$ (or in general in $C$) and increasing in the volatility aversion parameter $\gamma$.

Finally, regarding the influence of the factor model and the imposition of market neutrality in this new setting, note that $\mathrm{Rate}(t)$ is insensitive to it, since it only depends on the risk aversion parameter $\gamma$, the volatility of the residual process $\sigma\sigma^{\prime}$, and the transaction cost matrix $C$. Similarly, in $\mathrm{Aim}(t,x)$ it only appears through the term ${\mathbb{E}}_{t,x}[\mathrm{Frictionless}(s)]$ and hence only when considering the future optimal strategies in the frictionless case, which has been described in section 3. Likewise, the residual process $X_{s}$ only affects the strategy through this same term and hence, as seen in section 3 when studying these frictionless cases, only through the distance between this residual and its historical mean, like in classical pairs trading.

5 Monte Carlo simulations

We conclude the paper by providing some high-dimensional numerical simulations that give new insights about the behavior of the previously discussed strategies and their sensitivity to the different parameters, emphasizing in a separate simulation the role of the factor model and its connection with market-neutrality. To this end, we first simulate $X$ in the case of $N=100$ residuals by using exact Monte Carlo sampling along a discrete time grid, and we execute the previous strategies for some parametric choices of $X$ and some values of $p$ to compute sample paths of $\pi_{t}$ and $W_{t}$ and histograms of the terminal Profit & Loss (P&L). We have however opted to defer experiments with real data to a separate paper, since examining carefully the delicate issues of asset selection, rebalance frequency, construction of the factor models and obtention of $X$, high-dimensional parameter estimation and updating, risk control, etc. that the problem requires would be impossible to consider here without prohibitively extending the length of the paper.

During all this section, we therefore fix the following parameters for our model:

[TABLE]

$A$ is diagonal with entries drawn i.i.d from a normal distribution of mean 0.5 and standard deviation $0.1$, and the coefficients of $\sigma$ are drawn i.i.d from a uniform distribution in $[-0.3,0.3]$, except for the diagonal elements which are drawn from a uniform distribution in $[0,0.5]$. Furthermore $p=\mathds{1}$ for the first simulations, and we will also perturb it later to study different factor model regimes and the impact of imposing market neutrality. We also fix a finite horizon of $T=20$ and a temporal grid $0=t_{0}<t_{1}<\ldots t_{L}=T$ obtained by discretizing $[0,T]$ with constant $\Delta t=T/L=20/400=0.5$.

From a financial perspective, the above choice of parameters means that the 100 coordinates of $X$ are correlated and mean-revert with similar speeds (given by the eigenvalues of $A$) to an equilibrium of [math], describing an average number of approximately 10 oscillation cycles of ups and downs before the terminal time (given by the product of $T$ and the average mean-reversion speed). The choice of $p=\mathds{1}$ arises when the asset returns themselves are mean-reverting and may be modeled directly by $X$ so we can take $\Lambda=0$ in our factor model, while the perturbations of $p$ will imply departing from this assumption to factor models with heavier loadings, in which imposing market neutrality leads to more leveraged positions. As an illustration, we plot some sample paths of the first three coordinates of such a process $X$ in Figure 1 below.

We then sample $M=1,000$ paths of $X$ on this grid exactly with standard Monte Carlo techniques by using the fact that, since

[TABLE]

$X_{t+\Delta t}|X_{t}\sim N\left(\mu(X_{t},\Delta t),\Sigma(\Delta t)\right)$, where

[TABLE]

and execute the following strategies101010We have just simulated some simple cases of the previously discussed strategies for space limitation reasons, but it would be interesting as well to execute the strategies with some perturbations of the real parameters to simulate possible microstructural noise and imperfect estimation. at the corresponding times $t_{l}$’s, with $W_{0}=\pi_{0}=0$ and $\pi_{t}$ constant between consecutive times:

1. (1)

The exponential utility strategy of Theorem 3.1 with $\gamma(t)=1,2,3,4$ and $r=0,2\%$. 2. (2)

The mean-variance utility strategy with dollar-neutrality penalty of Corollary 4.1.1 with $\gamma(t)=1,2,3,4$, $\alpha(t)=0,20,50$ and $r=2\%$. 3. (3)

The mean-variance utility strategy with quadratic transaction costs of Theorem 4.1 with $\gamma(t)=1,2,3,4$, $r=2\%$, and $C=\lambda\sigma\sigma^{\prime}$ for $\lambda=0.1,0.5,1$.

Moreover, to study the result of imposing market-neutrality through the market-neutral portfolios constructed in section 2 under different factor model regimes, we perform the following additional simulation in which we experiment with the parameter $p$, which encapsulates all the factor model information and which we perturb to simulate the effect of going away from the case where the returns themselves are mean-reverting (which corresponds to the previous case $p=\mathds{1}$) and of having progressively more leveraged market-neutral portfolios:

1. (4)

The three strategies above with $\gamma(t)=1,\alpha(t)=0$ and $r=2\%$ (and $\lambda=1$ for the third strategy) for $p=\mathds{1}+\epsilon_{a}$ and $a=1,2,4,8$, where $\epsilon_{a}$ is a $N$-dimensional vector whose components are drawn i.i.d. from a uniform distribution in $[-a,a]$.

We present the simulated path of a sample wealth process $(W_{t})_{t\in[0,T]}$, the simulated path of the first coordinate of a sample allocation process $(\pi_{t})_{t\in[0,T]}$, and the histogram for the terminal wealth $W_{T}$ for each of the above cases in the following four subsections, along with a final analysis:

5.1 Simulations of the exponential-utility strategy

5.2 Simulations of the mean-variance strategy with dollar neutrality

5.3 Simulations of the mean-variance strategy with quadratic transaction costs

5.4 Simulations for different factor model and market-neutrality regimes

5.5 Comparison of the simulated strategies

We now present our main conclusions after observing the previous plots, analyzing the behavior of the histograms of the final wealth $W_{T}$, the sample paths of the wealth process $(W_{t})_{t\in[0,T]}$, and the sample paths of the positions $(\pi_{1t})_{t\in[0,T]}$, with a final subsubsection discussing the effect of imposing market-neutrality under different $p$’s.

5.5.1 Histograms of the final wealth

Looking first at the above histograms (Figures 3-4 (g)-(i), and 2 (e)-(f)), we see that, for our parametric choice and our setting in which $X$ is effectively a multidimensional Ornstein-Uhlenbeck process with known parameters,

1. (1)

The most profitable strategy is the one derived from the exponential utility (Figure 2, (e) and (f)) with the lowest risk-aversion parameter $\gamma$, even in the most adverse scenarios of the histogram and both with and without zero interest rates. Moreover, even for bigger values of $\gamma$ this strategy significantly performs better under any regime of $r$ and $\alpha$ than the mean-variance strategy (Figure 3, (g)-(i)). 2. (2)

We observe the following outcomes when changing one of the parameters for each of the three strategies (Figures 3-4 (g)-(i), and 2 (e)-(f)):

• •

Increasing the value of the risk-aversion parameter $\gamma$ produces a concentration of the density of $W_{T}$ around smaller values, i.e., the expected wealth decreases and so does the dispersion around it.

• •

Increasing the dollar-neutrality penalty $\alpha$ has this same negative effect, but makes little difference unless the increments in $\alpha$ are considerable.

• •

Increasing the value of the interest-rate $r$ has an overall positive effect, which is more pronounced in the mean-variance case since, as we mentioned at the end of section 2, the investor is then more aggressive than the exponential agent.

• •

Increasing the transaction cost parameter $\lambda$ decreases the expected terminal wealth, but it also skews its distribution producing a considerable right-tail (whereas all the other distributions are essentially symmetric).

These outcomes have a natural interpretation: since the model is perfectly specified and the parameters are known, the derived strategies will always produce benefits by construction, and they will be bigger the fewer additional constraints we impose (such as risk-aversion, dollar-neutrality, and transaction costs) and the more we can take advantage of previous success (by increasing $r$). This situation, however, might not apply under parameter misspecification, where the additional constrains would help the investor mitigate the model risk.

5.5.2 Sample paths of the wealth process $W_{t}$

Examining next the sample paths for the particular simulation which is plotted (Figures 3,4 (a)-(c), and 2 (a)-(b)), we observe the same patterns as discussed in the previous paragraph when modifying the parameters $\gamma,\alpha$, $r$ and $\gamma$. There are, however, two new observations:

1. (1)

In the three strategies, after an initial period of ups and downs and similarity between the different strategies, there is a tendency towards stabilization because of the asymptotic properties of the Ornstein-Uhlenbeck process, and of differentiation depending on the parametric choices. 2. (2)

This phenomenon is especially pronounced with the exponential utility and with bigger values of $r$ (Figure 2 (a)-(b)) since it takes more aggressive positions, reflecting the fact that sometimes the agent will invest more capital than what she will make at that moment (and sometimes even having temporary negative wealth and borrowing aggressively) to continue executing the strategy.

5.5.3 Sample paths of the positions $\pi_{1t}$

Considering now the plots of the sample paths of the positions $\pi_{1t}$ (Figures 3,4 (d)-(f), and 2 (c)-(d)), we similarly notice that

1. (1)

The positions become more extreme when decreasing $\gamma$, $\alpha$ and $\lambda$ (i.e., the risk-aversion parameter, the non-dollar-neutrality penalty and the transaction cost parameter) and when increasing $r$ (the interest rate). The greatest overall impact is produced by $\gamma$ and $\lambda$ and then $r$, especially in the mean-variance case for the same reasons as before. 2. (2)

The exponential utility strategy takes more extreme positions than the mean-variance strategies, which in this idealized setting of perfect estimation partially explains why the exponential agent obtains a greater wealth at the terminal time. 3. (3)

The cycles in the positions $\pi_{1t}$ match the oscillations of $X_{1t}$ depicted in Figure 1, as described theoretically in the corresponding equations.

5.5.4 Effects of imposing market neutrality

Finally, looking separately at the effect of imposing market neutrality under various factor model regimes depending on $p$ (which, as we mentioned, depends quadratically on the factor model loadings), we observe the following (Figure 5):

1. (1)

As the parameter $p$ gets bigger, the market neutral portfolios of section 2 become more extreme and the adopted positions $\pi_{t}$ also become more aggressive, especially in the exponential utility case (Figure 5, (d)-(f)). 2. (2)

Since the strategy is more aggressive but we have perfect estimation, with bigger $p$ the mean-variance and especially the exponential strategy become more profitable. However, the wealth process also has more ups and downs (Figure 5, (a)-(c)), the standard deviation of the terminal wealth increases considerably (Figure 5, (g)-(h)), especially in the mean-variance case, and with the biggest $p$ there are also heavy losses when transaction costs are incorporated (Figure 5, (c),(i)). The strategies are therefore riskier, but a relatively large value of $p$ is needed to appreciate its effect. 3. (3)

Lastly, note that the influence of $p$ on the strategies also depends most of the time on the value of $r$, since they normally appear combined as a factor of $rp$ in the equations describing the strategies. In particular, when $r=0$ there is no theoretical effect associated to $p$ (apart from possible model risk and high leverage in a real-world setting) unless the dollar-neutrality parameter $\alpha(t)\neq 0$.

6 Conclusions and further research

In this paper we have aimed to provide a systematic study of high-dimensional statistical arbitrage combining stochastic control and factor models. To this end, we have first proposed a general framework based on a statistically-constructed factor model, and then shown how to obtain analytically explicit market-neutral portfolios and rephrase our problem in terms of them to make it tractable and get market neutrality by construction. Using this insight, we have studied the question of optimizing the expected utility of the investor’s terminal wealth in continuous time under both an exponential and a mean-variance objective. In both cases, we have obtained explicit closed-form solutions that avoid potentially unfeasible high-dimensional numerical methods, analyzed the corresponding strategies from the perspective of statistical arbitrage and the underlying factor model, and discussed extensions involving the addition of soft constraints on the admissible portfolios (like dollar-neutrality) and the presence of temporary quadratic transaction costs. Finally, we have run Monte Carlo simulations to explore the behavior of the previous strategies, and analyzed their qualitative aspects and their sensitivity to the relevant parameters and the underlying factor model.

There are four natural extensions to our work, on which we are conducting research at the moment and which we intend to publish in separate papers. First, one could investigate a more realistic version of the problem in which, rather than in continuous time, the investor may only trade more realistically at an increasing sequence of optimally chosen stopping times,generalizing the literature initiated by Leung and Li (2015) and developing robust and efficient numerical methods. Second, it would be interesting to study more realistic modelizations of market frictions, illiquidity, and transaction costs, or to develop a model considering issues of parameter misspecification. Third, on a more empirical side and as we mentioned at the start of the section 5, one should consider in this setting the problems of construction of the factor models, high-dimensional parameter estimation, and risk control, along with out-of-sample experiments with real market data under the strategies developed in this paper. Fourth and finally, one could study a more data-driven version of the problem, where the fixed stochastic model is replaced by new tools from reinforcement learning.

7 Acknowledgements

The author would like to thank George Papanicolaou for suggesting the topic of the previous research and for many insightful discussions about the problem and the presentation of the results, and the editor and an anonymous reviewer for their very helpful suggestions to improve the quality of the paper.

Appendix A Appendix: Proofs

A.1 Proof of Theorem 3.1

**Proof of Proposition 3.1:

**The dynamic programming principle suggests that the value function $H$ should satisfy the dynamic programming equation (2.5) with terminal condition $H(T,x,w)=-e^{-\gamma w}$. The optimal control may then be found in feedback form by looking at the first order condition of the term inside the supremum, since the corresponding function is quadratic and concave in $\pi$ (if $\partial_{ww}H<0$, i.e., if there is risk aversion). The first order condition gives that

[TABLE]

and solving for $\pi$ we find the control given in the proposition’s statement. Putting it back into (2.5) we get the following non-linear and $(N+2)$-dimensional PDE

[TABLE]

Now, looking at the terminal condition, we guess that the solution of this PDE will be of the form $H(t,x,w)=-\exp(-\gamma(we^{r(T-t)}+h(t,x)))$ for some function $h(t,x)$ to be determined and such that $h(T,x)=0$. Some easy computations then show that

[TABLE]

[TABLE]

Plugging all this into (A.1), dividing everything by $-\gamma H$, and doing some simple algebra to expand the last term yields

[TABLE]

and we can see that the non-linear terms in $h$, the terms in $w$, and the third and part of the last term of the PDE get cancelled and the equation gets considerably simplified, obtaining

[TABLE]

This is now a parabolic linear PDE in $h$ and we can find explicitly its solution by using the Feynman-Kac formula. Indeed, if we consider the stochastic process given by

[TABLE]

we can rewrite the above equation in terms of the infinitesimal generator $\mathcal{L}^{*}$ of $Y$ as

[TABLE]

and then we can express its solution via the following conditional expectation, which is the probabilistic representation given in the proposition’s statement:

[TABLE]

Finally, to find $h$ explicitly, notice that we can easily solve the SDE (A.2), obtaining, for $s\geq t$,

[TABLE]

and this allows us to compute the two expectations in our expression for $h$ above. Indeed, Fubini’s theorem and elementary facts about the Brownian motion immediately yield

[TABLE]

and, interchanging integral and expectation again and noticing that

[TABLE]

we similarly find out that

[TABLE]

which gives us the complete explicit solution of the DPE, and hence the explicit form of the optimal strategy $\pi^{*}$ by using equation (3.1). ∎

**Proof of Proposition 3.2:

**Since in the previous proof we have found explicitly the classical smooth solution $H$ of the dynamic programming equation, we just have to check that $\pi^{*}\in\mathcal{A}_{[0,T]}$ and that the usual regularity conditions hold for the classical proof to apply. More precisely, this means that the local martingale $dH-\mathcal{L}_{t,x,w}^{\pi}Hdt$ is a supermartingale for any admissible $\pi$ and a true martingale for $\pi^{*}$, where $\mathcal{L}_{t,x,w}^{\pi}$ is the infinitesimal generator of the controlled process $(X,W^{\pi})$, or some sufficient condition for this like the one we stated in Proposition 3.2 in terms of the model parameters, which is what we will show here.

As for the first issue, it is easy to see that $\pi^{*}\in\mathcal{A}_{[0,T]}$. Indeed, it is obviously $\mathcal{F}_{t}$-adapted and predictable (in fact, it has continuous paths) and, using the trivial inequalities $||x+y||^{2}\leq 2||x||^{2}+2||y||^{2}$ and $||Ax||\leq||A||||x||$ and the fact that $X_{t}$ is a Gaussian process, it is easy to see that $\int_{0}^{T}{\mathbb{E}}[||\pi_{s}^{*}||^{2}]ds<\infty$. Moreover, applying Ito’s formula to the process $e^{-rt}W_{t}$ yields

[TABLE]

and, therefore,

[TABLE]

for any $t\geq 0$ and any admissible control $\pi$. Thus, the SDE for $W$ has a unique strong solution $W^{*}$ for the particular case $\pi=\pi^{*}$ for any initial data, given by the above integral (note that the stochastic integral is well defined, since $dX_{s}=A(\mu-X_{s})ds+\sigma dB_{s}$, $\pi^{*}$ and $X$ are continuous, and again $\int_{0}^{T}e^{-2rs}{\mathbb{E}}[||\pi_{s}^{\prime*}\sigma||^{2}]ds<\infty$).

As for the regularity conditions, we adapt the proof of Theorem 2.1. of Lintilhac and Tourin (2016), which guarantee the uniform ${\mathbb{P}}$-integrability of the family of random variables $(H(\tau,X_{\tau},W_{\tau}^{*}))_{\tau\in[0,T]}$ where $\tau$ is a $\mathcal{F}$-stopping time, and which we adapt to the parameters of the present model obtaining the sufficient conditions stated in Proposition 3.2.

The key observation to adapt their proof is that in our case we also have that the hypothesized value function is of the form $H(t,x,w)=-\exp(-\gamma we^{r(T-t)}-\frac{1}{2}x^{\prime}A_{2}(t)x-A_{1}(t)x-A_{0}(t))$ for some explicit smooth functions $A_{i}(t)$ that we computed in the proof of Proposition 3.1, and our $X$ is also a matrix Ornstein-Uhlenbeck process under ${\mathbb{P}}$ with SDE $dX_{t}=A(\mu-X_{t})dt+\sigma dB_{t}$, and

[TABLE]

as we showed in (A.3). Thus, using the Cauchy-Schwarz inequality as in their proof, the part corresponding to $-\frac{1}{2}X_{\tau}^{\prime}A_{2}(\tau)X_{\tau}-A_{1}(\tau)X_{\tau}-A_{0}(\tau)$ in the above expression for $H(\tau,X_{\tau},W^{*}_{\tau})$ may be bounded following their reasoning. As for the part corresponding to $-\gamma W_{\tau}^{*}e^{r(T-\tau)}$, we can again repeat their argument, but noting that the quadratic term in $X_{s}$ in the first integral above is now $X_{s}^{\prime}C_{0}(s)X_{s}$ for the matrix $C_{0}(s)$ that we defined before, and likewise the term in $X_{s}$ in the second integral is $X_{s}^{\prime}C_{1}(s)$, which following their proof gives respectively the two explicit sufficient conditions that we stated in Proposition 3.2. ∎

A.2 Proof of Theorem 3.2

The proof of this follows the same lines as the previous one and is simpler, so we just indicate the relevant changes. The HJB equation is now

[TABLE]

with terminal condition $H(T,x,w)=w.$

Guessing that the value function will now be of the form $H(t,x,w)=we^{r(T-t)}+a(t)+b(t)^{\prime}x+\frac{1}{2}x^{\prime}c(t)x$ for a scalar $a(t)$, an $N$-dimensional vector $b(t)$, and a symmetric $N\times N$ matrix $c(t)$, and plugging this into the above equation, we obtain the hypothesized optimal control given in the statement of the theorem and the above PDE gets reduced to the following system of three first-order linear matrix ODEs

[TABLE]

[TABLE]

[TABLE]

with terminal conditions $a(T)=b(T)=c(T)=0$.

This system has an explicit bounded solution in $[0,T]$, since the classical solution of the general first-order linear matrix ODE $\partial_{t}y+uy+v(t)=0$ with $y(T)=0$ is given by

[TABLE]

if $v(s)$ is continuous on $[0,T]$, in which case it is automatically bounded on $[0,T]$ as well; and similarly the classical solution of $\partial_{t}y+uy+yu^{\prime}+v(t)=0$ with $y(T)=0$ for a symmetric $v$ is given by

[TABLE]

Thus, the HJB equation has an explicit classical solution which has quadratic growth in the state variables uniformly in $t$. A classical verification result (cf. for example Theorem 4.3 of Guyon and Labordère (2013)) then yields that our hypothesized optimal control is indeed optimal provided that it is admissible, which may be checked exactly as in the proof of Theorem 3.1. ∎

A.3 Proof of Theorem 4.1

**Proof of Proposition 4.1:

**The corresponding dynamic programming equation is in this case

[TABLE]

with terminal condition $H(T,x,w,\pi)=w$ and where the supremum is obviously attained at $I^{*}=C^{-1}\nabla_{\pi}H$.

Substituting this back in the above equation and plugging the stated ansatz we obtain that

[TABLE]

Matching the coefficients for $\pi^{\prime}(\cdot)\pi$, $\pi^{\prime}(\cdot)$, and the constant yields the above differential equations. ∎

Before we prove the next proposition, we state here the following result for comparison and existence of solutions of matrix Riccati ODEs (cf. Theorem 2.2.2 in Kratz (2011)), which we will use in our proof.

Theorem A.1**.**

Let $L_{1}(t),L_{2}(t),L(t),N_{1}(t),N_{2}(t)\in{\mathbb{R}}^{d\times d}$ be piecewise continuous on ${\mathbb{R}}$. Moreover, suppose $L_{1}(t),L_{2}(t),N_{1}(t),N_{2}(t)$ and $S_{1},S_{2}\in{\mathbb{R}}^{d\times d}$ are symmetric. Let $T>0$ and

[TABLE]

on $[0,T]$. Assume that the terminal value problem

[TABLE]

has a (symmetric) solution $H_{1}$ on $[0,T]$. Then the terminal value problem

[TABLE]

has a (symmetric) solution $H_{2}$ on $[0,T]$ and $H_{1}(t)\geq H_{2}(t)$ for all $t\in[0,T].$

We are now in a position to give the following:

Proof of Proposition 4.2:

1. (1)

The first statement follows directly from the comparison Theorem A.1 stated before, since the matrix Riccati ODE

[TABLE]

with terminal condition $a(T)=0$ has the obvious symmetric solution $a(t)=0$ defined on all $[0,T]$. Thus, (4.2) has a symmetric classical solution $a(t)$ defined on all $[0,T]$ with $a\leq 0$, which is bounded because $[0,T]$ is compact and $a$ is differentiable hence continuous.

As for the particular solution when $\gamma(t)=\gamma$, simply note that pre- and post-multiplying (4.2) by $C^{-1/2}$ and defining $\tilde{a}:=C^{-1/2}aC^{-1/2}$ gives the new Riccati

[TABLE]

whose solution is $\tilde{a}(t)=D\tanh(D(t-T))$. 2. (2)

The existence of solutions with polynomial growth and their probabilistic representation in the above form follow from a vector-valued version of the Feynman-Kac theorem (see Appendix A.3 of Cartea et al. (2018) for a proof of how to adapt the one-dimensional case) provided that the appropriate regularity conditions hold. Using, for example, Condition 2 of Appendix E in Duffie (2010), it is sufficient that all the functions of $(t,x)$ $A(\mu-x)$, $\sigma$, $a(t)^{\prime}C^{-1}$, $e^{r(T-t)}(A(\mu-x)-rp)$ (and $b(t,x)^{\prime}C^{-1}b(t,x)$ for the existence of $d$) are uniformly Lipschitz in $x$, they and their first and second derivatives in $x$ are continuous with polynomial growth in $x$ uniformly in $t$, and $a(t)\leq 0$. All of these properties are straightforward to check in this case because all the corresponding functions are given explicitly and are simple, and the required properties for $a$ follow from (1).

The fact that $b$ has linear growth in $x$ uniformly in $t$ is then a consequence of the probabilistic representation (4.5). Indeed, Fubini’s theorem implies that

[TABLE]

whereas the fact that

[TABLE]

yields

[TABLE]

Combining the two pieces and using the boundedness of $a$ and the compactness of $[0,T]$ gives the desired uniform bound in $t$.

The quadratic growth of $d$ in $x$ uniformly in $t$ is then obvious looking at its probabilistic representation and using the linear growth of $b$.∎

**Proof of Proposition 4.3:

**Combining the two previous propositions, we have already found an explicit classical solution of the associated HJB equation with quadratic growth in the state variables uniformly in $t$, so using again Theorem 4.3 in Guyon and Labordère (2013), we just have to verify that the candidate intensity given in Proposition 4.1 is admissible.

For this, first of all note that the corresponding SDEs controlled by the above intensity have a unique strong solution for any initial data. Indeed, given $I^{*}$ and the definition of $I$ as $d\pi=Idt$, we can solve explicitly the corresponding first-order linear matrix ODE for $\pi^{*}$ yielding, for $s\geq t$,

[TABLE]

and this $\pi^{*}$ in turn defines $W^{*}$ like in the proof of Theorem 3.1.

Finally, from the above construction it is obvious that both $\pi^{*}_{t}$ and $I^{*}_{t}$ are $\mathcal{F}_{t}$-adapted and predictable (in fact, they have continuous paths), and the property that $\pi^{*}$ is in $L^{2}([0,T]\times\Omega)$ (i.e., that $\int_{0}^{T}{\mathbb{E}}[||\pi_{s}^{*}||^{2}]ds<\infty$) stems from the observation that $\mathrm{Rate}(u)$ is deterministic and bounded (because of Proposition 4.2.1), $\mathrm{Aim}(t,x)$ has linear growth in $x$ uniformly in $t$ (by Proposition 4.2.2), and $X$ is a Gaussian process (so it is in $L^{2}([0,T]\times\Omega)$).

$I_{t}^{*}$ is likewise in $L^{2}([0,T]\times\Omega)$ since, as we saw in Proposition 4.1, $I_{t}^{*}=C^{-1}(a(t)\pi_{t}^{*}+b(t,X_{t}))$ and we can therefore use the triangular inequality, the just shown fact that $\pi_{t}^{*}$ is in $L^{2}([0,T]\times\Omega)$, and the same arguments as above that $a(t)$ is bounded (because of Proposition 4.2.1), that $b(t,x)$ has linear growth in $x$ uniformly in $t$ (by Proposition 4.2.2), and that $X$ is a Gaussian process, to conclude. ∎