Portfolio Optimization for Cointelated Pairs: SDEs vs. Machine Learning
Babak Mahdavi-Damghani Konul Mustafayeva Stephen Roberts Cristin Buescu





 
 a4paper,
 
 left=28mm,
 right=28mm,
 top=30mm,
 bottom=28mm
 









 

scopes,matrix,positioning















 

⌈⌉
⌊⌋







      
theoremTheorem
assumptionAssumption
lemmaLemma

remarkRemark
propositionProposition
conjectureConjecture 
questionQuestion 
definitionDefinition
exampleExample
corollaryCorollary





iblabel[1]





  mymx/.style=matrix of math nodes,nodes=myball,column sep=4em,row sep=-1ex,
  myball/.style=draw,circle,inner sep=0pt,
  mylabel/.style=midway,sloped,fill=white,inner sep=1pt,outer sep=1pt,below,
    execute at begin node=,execute at end node=,
  plain/.style=draw=none,fill=none,
  sel/.append style=fill=green!10,
  prevsel/.append style=fill=red!10,
  route/.style=-latex,thick,
  selroute/.style=route,blue!50!green   cirwhite/.style=draw=gray,circle,fill=white,minimum size=1pt,inner sep=2pt,line width=0.2mm,
            cirred/.style=draw=gray,circle,fill=red,minimum size=1pt,inner sep=2pt,line width=0.2mm,

1]Babak Mahdavi-Damghani2]Konul Mustafayeva2]Cristin Buescu1]Stephen Roberts[1]Oxford-Man Institute of Quantitative Finance, Oxford, UK[2]Department of Mathematics, King's College London, London, UK 


Portfolio Optimization for Cointelated Pairs: SDEs vs Machine Learning

    [
    

=======================================================================

  
With the recent rise of Machine Learning as a candidate to partially replace classic Financial Mathematics methodologies, we investigate the performances of both in solving the problem of dynamic portfolio optimization in continuous-time, finite-horizon setting for a portfolio of two assets that are intertwined. 

In Financial Mathematics approach we model the asset prices not via the common approaches used in pairs trading such as a high correlation or cointegration, but with the cointelation model in <cit.> that aims to reconcile both short-term risk and long-term equilibrium. We maximize the overall P&L with Financial Mathematics approach that dynamically switches between a mean-variance optimal strategy and a power utility maximizing strategy. We use a stochastic control formulation of the problem of power utility maximization and solve numerically the resulting HJB equation with the Deep Galerkin method introduced in <cit.>. 

We turn to Machine Learning for the same P&L maximization problem and use clustering analysis to devise bands, combined with in-band optimization. Although this approach is model agnostic, results obtained with data simulated from the same cointelation model as FM give an edge to ML. 

 Correspondence to: Konul Mustafayeva, Department of Mathematics, King's College London, Strand, London,  WC2R 2LS. 

E-mail: konul.mustafayeva@kcl.ac.uk 

81

empty

§ INTRODUCTION
 
In a financial market with two assets that exhibit clear dependence we go beyond high correlation which is used in pairs trading, and model asset prices with the hybrid model of cointelation introduced in <cit.>. We solve the portfolio optimization problem employing a more general set of admissible strategies than long/short strategies used in pairs trading. 

A pairs trading strategy involves matching a long position with a short position in two assets with a high correlation. Pairs trading was pioneered in the mid 1980s by a group of quantitative researchers from Morgan Stanley (for introduction to pairs trading see <cit.>). The securities in a pairs trade must have a high positive correlation, which is the primary driver behind the strategy’s profits. 

Pairs trading is based on the high historical correlation of two assets and a trader's view that the two securities will maintain a specified correlation. A pairs trading strategy is applied when a trader identifies a correlation discrepancy. More specifically, the trader monitors performance of two historically correlated securities. When the correlation between the two securities temporarily weakens, i.e. the spread widens, the trader applies a trading strategy which shorts the high asset and buys the low asset. As the spread narrows again to some equilibrium value, a profit results. 

However, many authors argue that correlation is an inappropriate measure of dependency in financial markets, since returns often exhibit a nonlinear co-dependence (e.g. <cit.>, <cit.>). Mahdavi-Damghani, et al. <cit.> showed that measured correlation of the returns of a mean-reverting processes is misleading: a strong positive correlation does not necessarily imply that two stochastic processes move in the same direction and vice versa. Cointegration, on the other hand, tests the long-term equilibrium relationships between assets and has been extensively used in pairs trading ( see <cit.>). Cointegration tests do not measure how well two variables move together, but rather whether the difference between their means remains constant. Sometimes series with high correlation will also be cointegrated, and vice versa, but this is not always the case.

The cointelation model was introduced in <cit.> as a hybrid model which reconciles correlation and cointegration by capturing both short-term risk and long-term equilibrium. The rationale for the long term risk is that during the time of rare market crashes all assets prices fall. However, in the more bullish periods, the short term risk increases, the long term risk becomes less pronounced and the “macro” driver less visible. These influences are accompanied with mean reversion forces from one asset to the other.  

In this setting we consider a continuous-time, finite horizon portfolio optimization problems for pairs of assets whose prices follow the cointelation model in <cit.>. Generally, the optimization problem is to find the optimal control

    w̃^* = w̃∈𝒜argmax  U(X^w̃_t, Y^w̃_t)

where U(x) is a utility function, w̃ = (w̃_1, w̃_2) is a vector of proportions of wealth invested in each asset, 𝒜 is a set of admissible strategies: either w̃_1 = - w̃_2 (long/short) or w̃_1, w̃_2 > 0 with w̃_1 + w̃_2 = 1 (long only).

We solve the portfolio optimization problem in (<ref>) with Financial Mathematics and Machine Learning methodologies and compare their performance. In Financial Mathematics approach we use SDE evolution of asset prices, whereas the Machine Learning approach does not assume an underlying model and applies generally to any pair of assets. 

In Section <ref> we review the cointelation model. In Section <ref> we use the classical Financial Mathematics criteria: mean-variance optimization and power utility maximization. In Section <ref> we use clustering analysis from Machine Learning to solve the P&L maximization problem. We present the results of each approach in Section <ref> and discuss them comparatively.
 


§ REVIEW OF COINTELATION MODEL FOR PAIRS OF ASSET 

We first present the usual way correlation is calculated in the financial industry (see e.g. p.274 <cit.>, <cit.>). Assume we have two assets with prices modeled by stochastic processes (X_t)_t≥ 0 and (Y_t)_t≥ 0 on a probability space (Ω , ℱ, ℙ). We have N observations of X and Y at intervals Δ t, i.e. X(t_i) and Y(t_i) with i=1,...,N and Δ t = t_i - t_i-1. Here  Δ t ∈{1, 5, 22, 252} corresponds to daily, weekly, monthly and yearly data. The Δ t-returns on i-th data point of assets X and Y is

    R_X(t_i, Δ t)    =   X(t_i + Δ t)  - X(t_i)/X(t_i)
    
        R_Y(t_i, Δ t)    =   Y(t_i + Δ t) - Y(t_i)/Y(t_i).

The sample volatilities of time series of asset prices X and Y are then

    σ_X(Δ t)    =   √(1/Δ t(N-1)∑_i=1^N (R_X(t_i, Δ t) - R_X)^2)
    σ_Y(Δ t)    =   √(1/Δ t(N-1)∑_i=1^N (R_Y(t_i, Δ t) - R_Y)^2),

where R_X, R_Y are the sample average of all the returns in the series of X  and Y, respectively. 

The sample covariance between the returns of assets X and Y is given by

    σ_XY(Δ t) = 1/Δ t(N-1)∑_i=1^N (R_X(t_i, Δ t) - R_X)(R_Y(t_i, Δ t) - R_Y).

In this paper we consider the measured correlation, which is the sample cross-correlation given by

    ρ_XY(Δ t) = σ_XY(Δ t)/σ_X(Δ t) σ_Y(Δ t).


For correlation to be an appropriate choice of measure of co-dependence the assumption of linear dependency between series needs to be satisfied (see Chapter 1.4 <cit.>). Often in financial markets with a non-linear dependence between returns, the correlation is an inappropriate measure of co-dependency and is misleading, especially when used to capture long-term relationship between assets (see <cit.> and <cit.>).  

An alternative statistical measure to correlation is cointegration. If two time series X_t and Y_t are integrated[A time series X_t is integrated of order d if (1-L)^dX_t is a stationary process. Here L is a lag operator.] of order d and there exists β such that a linear combination X_t + β Y_t is integrated of order less that d, then X_t and Y_t are cointegrated (see <cit.>). Since the spread of cointegrated asset prices is mean reverting, they have a common stochastic trend, i.e. the asset prices are ‘tied together’ in the long term, although they might drift apart in the short-term (see <cit.>). Because the cointegration requires sophisticated statistical analysis, it has not been used as widely as correlation in the financial industry. 






Although correlation and cointegration are related, they are different concepts. High correlation does not necessarily imply high cointegration, and neither does high cointegration imply high correlation (e.g. see Figure 4 in <cit.>). Two assets may be perfectly correlated over short timescales yet diverge in the long run, with one growing and the other decaying. Conversely, two assets may follow each other, with a certain finite spread, but with any correlation, positive, negative or varying. 

Mahdavi-Damghani <cit.> proposed cointelation as a hybrid model that aims to mediate between correlation and cointegration. It captures both short-term and long-terms relationships between the assets. 

 Consider a filtered probability space by (Ω, ℱ, (ℱ_t)_(t ≥ 0), ℙ), with the historical probability measure, ℙ. The cointelation model for a pairs of assets with prices X_t and Y_t defined in <cit.> as 

    dX_t = μ X_t dt + σ X_t dW_t,   
           dY_t = κ(X_t-Y_t)dt  + η Y_t dW̃_t, 
           d ⟨ W, W̃⟩_t = ρ dt,

where μ∈ℝ, σ>0, X(t_0) = x_0 are the drift, diffusion coefficients and initial value of asset price X; 0< κ≤ 1, η > 0, Y(t_0) = y_0 > 0 are the rate of mean reversion, volatility and initial value of the asset price Y; (W̃(t))_t≥ 0 and (W(t))_t≥ 0 are two correlated Brownian motions with constant correlation coefficient -1 ≤ρ≤ 1 that generate the filtration (ℱ_t)_(t ≥ 0).

The processes (X)_t≥ 0 and (Y)_t≥ 0 are called the leading process and the lagging process, respectively. This is due to the fact that the lagging process reverts around the leading process.

We present here the concepts of inferred correlation function and number of crosses formula introduced in <cit.> in order to device a test whether two pairs are cointelated .  

Let ρ_XY^*(Δ t) be the inferred correlation function between two times series of cointelated asset prices defined as follows

    ρ^*_XY(Δ t) = sup_0<Δ̃t≤Δ tρ_XY(Δ̃ t).

Sometimes there may not be enough data to calculate Δ t-inferred (measured) correlation of cointelated assets. In <cit.> the following formula for approximation of inferred correlation (<ref>) was proposed via examining various data sets:

    ρ^*_XY(Δ t) ≈ρ + (1-ρ) [1 - exp(- λκ(Δ t -1))],

where  κ∈ [0,1], λ >0, ρ∈ [-1,1]. The parameter λ≈ 1.75 for "regular financial data", although it is itself a function in general. Thus, if one does not have enough empirical data to calculate, for example, the yearly (252 days) inferred correlation, the formula in equation (<ref>) allows to approximate it using only κ and ρ parameters of cointelation model in (<ref>) and setting Δ t =252, λ≈ 1.75.

The motivation for inferred correlation approximation formula (<ref>), is that in the discrete version of the processes in equation (<ref>) the measured correlation increases as the time increment, Δ t, increases (e.g. correlations calculated using daily, weekly, monthly returns). Moreover, the measured correlation of cointelated pairs will converge to 1 faster as the speed of mean reversion parameter κ increases. If we set ρ=-1 in (<ref>), the inferred correlation of cointelated asset prices may cover the whole correlation spectrum [-1,1] (see Figure <ref>).

Another way for testing if two times series are cointelated is to study how many times the normalized series cross paths. If one discretizes equation (<ref>), then one can approximate the expectation of the number of times, Γ_x,y, the two stochastic process, x=X_i∈[1,2,… N] and y=Y_i∈[1,2,… N], cross paths as follows

    E[Γ_x,y(κ,N)] ≈  N[γ(1-κ) + 1/2√(κ)]

with N is the length of the data, γ is a positive constant and κ is the speed of mean reversion in equation (<ref>).

Compared to the number of times purely correlated SDEs (eg: without the mean reversion component, i.e. when κ = 0) the number of times the discrete version of the cointelated SDEs cross paths is larger than if they were random, and the bigger the κ the more often the paths of discretized SDEs cross each other per unit of time. 

Then two stochastic processes are cointelated (see <cit.>) if

  *  Inferred correlation formula in equation (<ref>) is verified;
    
  *  The number of crosses formula in equation (<ref>) is verified;
    
  *  the underlying assets have a reasonable physical connection that would suggest their spread should mean revert (e.g. oil and BP share prices). 


The parameters in cointelation model (<ref>) can be estimated using the inferred correlation formula (<ref>) and the number of crosses formula (<ref>) (see <cit.>). Similarly to the variance reduction methodology described in <cit.>, <cit.>, we can define 

    B_+   =   |max(X_t-Y_t,t∈[0,T])/2|, 
    
    B_-   =   |inf(X_t- Y_t,t∈[0,T])/2|.

We note that the estimation of κ has a higher variance when 
    Z_ρ = B_+>|X_t-Y_t|>B_-,

where ρ, on the other hand has quality samples. The reverse is true when

    Z_κ=|X_t-Y_t|>B_+⋃ |X_t-Y_t|<B_-.

We can therefore sample  κ in Z_κ and ρ in Z_ρ. Figure <ref> illustrates this.















 



§ FINANCIAL MATHEMATICS APPROACH FOR PORTFOLIO OPTIMIZATION PROBLEM

We consider the portfolio of two assets and model the their prices with the cointelation in (<ref>). We approach the optimization problem of this portfolio with classic Financial Mathematics criteria: mean-variance and power utility maximization. Since the cointelated assets are characterized by both correlation and mean-reversion components, we formulate the mean-variance optimization problem for long only strategies and we calculate the optimal strategies to make profit on correlation. To make profit on mean-reversion property of the cointelated assets we use stochastic control formulation of the power utility maximization problem for long/short strategies and calculate the optimal weights. We then maximize portfolio P&L by dynamically switching between these two optimal strategies.  




 §.§ Mean-variance optimization

We first review fundamental notions and concepts for mean-variance optimization.


  
Returns:
A portfolio considers a combination of n potential assets, with an initial capital V(0) and weights w_1, w_2,..., w_n, such that ∑_i^n w_i =1, w_iV(0) is the amount invested in security i for i= 1,2,...,n at time t=0.  
The number of shares to invest in security i at time t=0 is 

    n_i = w_iV(0)/S_i(0).

The value of portfolio at time t is

    V(t) = ∑_i=1^N n_i S_i(t).

Given the number of shares n_i with i = 1,...,n, the percentage of the portfolio invested in asset i at time t is

    w_i(t) = n_iS_i(t)/∑_i=1^N n_i S_i(t),

with ∑_i=1^N w_i(t) = 1.
The rate of return of asset i at time t (i.e. over [t -Δ t,t]) is given by

    R_i(t) = S_i(t) - S_i(t-Δ t)/S_i(t-Δ t) = S_i(t)/S_i(t-Δ t) - 1.

The rate of return of portfolio, R_p(t), is then

    R_p(t) = V(t) - V(t-Δ t)/V(t-Δ t).

We can show that the return of portfolio is a linear combination of the returns of individual assets as follows

    R_p(t)    =   V(t)/V(t-Δ t) - 1 = ∑_i=1^N n_i S_i(t)/∑_j=1^N n_i S_i(t-Δ t) - 1 
       =   ∑_i=1^N n_i S_i(t-Δ t)S_i(t)/∑_j=1^N n_i S_i(t-Δ t)S_i(t-Δ t) - 1 = ∑_i=1^N w_i(t)(R_i(t)+1) - 1 
       =   ∑_i=1^N w_i(t)R_i(t).

Sometimes it is more convenient to use log returns, which are defined for asset i by 

    r_i(t) = ln(S_i(t)/S_i(t-Δ t)).

It should be pointed out that for short period of time the log return is approximately equal to the rate of return 

    r_i(t) = ln(S_i(t)/S_i(t-Δ t)) = ln(R_i(t) +1) ≈ R_i(t).

Therefore we do not distinguish between these two returns, as long as the time increment, Δ t, is short. Going forward we will use daily logarithmic returns. Thus, the return of portfolio, r_p, at time at time t in this case becomes

    r_p = ∑_i=1^N w_ir_i.


  
Expectation and variance of returns:
By the linearity property of expected value operator, the expected return of portfolio, E(r_p), is 

    E(r_p) = E( ∑_i=1^N w_i r_i ) = ∑_i=1^N w_i E(r_i) = ∑_i=1^N w_i μ_i = w^⊤μ ,

where μ_i denotes the expected return of asset i and w^⊤ = [w_1, w_2, ..., w_n], μ = [μ_1, μ_2, ..., μ_n]^⊤.

The variance of the return of portfolio, Var(r_p), is given by

    Var(r_p)    =    E [ (∑_i=1^N w_i r_i - E(r_p))^2  ] = E [( ∑_i=1^N w_i(r_i - E(r_i)) )^2 ]
       =    E [ ( ∑_i=1^N w_i(r_i - E(r_i)) ) (  ∑_j=1^N w_j(r_j - E(r_j) )) ] 
       =   ∑_i=1^N ∑_j=1^N w_iw_j E[ (r_i - E(r_i))(r_j - E(r_j))]_:= σ(r_i,r_j)
       =   ∑_i=1^N ∑_j=1^N w_iw_j σ(r_i,r_j) = w^⊤Σ w,
 
where Σ denotes the covariance matrix of the asset returns, composed of all covariances between the returns of assets i and j defined as  σ(r_i,r_j). The variance of asset i’s return, which 
constitute the diagonal of the covariance matrix, is σ(r_i,r_i).



  §.§.§ Optimal investment strategy using mean-variance criterion

We consider a portfolio consisting of two assets. The uncertainty is modelled by a probability space (Ω, ℱ, P) with a filtration (ℱ_t)_t≥ 0 generated by two-dimensional Brownian
motion: (W, W̃). Denote by X(t) and Y(t) the prices of two assets at time t, with dynamics following cointelation model in (<ref>).
The investment behavior is modelled by an investment strategy h = (h_1, h_2). Here, h_i ∈ [0,1], i = 1, 2, denotes the percentage of total wealth invested in i-th asset (see equation (<ref>)). Let h_1(t) and h_2(t) denote respectively the portfolio weights for assets X and Y at time t. The holdings are allowed to be adjusted continuously up to a fixed horizon T.  

Denoting by V^h_t the value of portfolio at time t associated to a strategy h we have

    V^h(t) = h_1(t)V^h(t)/X(t)X(t) + h_2(t)V^h(t)/Y(t)Y(t) ,

with initial wealth V^h(t_0) = v_0. We restrict our considerations to self-financing strategies, where the value of the portfolio changes only because the asset prices change, i.e. there is no inflow or withdrawal of money <cit.>. In this case the dynamic of the wealth process is 

    dV^h(t)=V^h(t)[ h_1(t)dX(t)/X(t) + h_2(t)dY(t)/Y(t)].


Let 𝒜^1 denote the set of all admissible strategies, h = (h_1, h_2), satisfying:

  (i) Given v_0 > 0 the wealth process V^v_0,h(·) corresponding to w_0,h satisfies 

    V^v_0,h(t) ≥ 0,    0 ≤ t ≤ T,

  (ii)h_i(t) ≥ 0 for all i=1,2,

  (iii)∑_i=1^2h_i(t) = 1. 


An investment strategy, h ∈𝒜^1, is called optimal if there exists no other strategy h̃∈𝒜^1 such that E(r_p(h)) ≥ E(r_p(h)) and Var(r(h)) ≤ Var(r(h)) with at least one inequality being strict (see <cit.>).

We define a utility function, U(t,h), as in <cit.>:

    U(t,h) = 2τ E[r_p(t)] - σ^2[r_p(t)],

where τ≥ 0 is the risk tolerance coefficient. Then according to <cit.> we have the following proposition.

Finding an optimal strategy for mean-variance criteria is equivalent to the utility maximization problem:

    max_h(t) U(t, h)

with constraints

  * ∑_i=1^N h_i = 1,

  * h_i ≥ 0   ∀ i.

and U(t, h) given in (<ref>).

Thus we have optimization problem in equation (<ref>). 
From equation (<ref>) we have that the rate of return of our portfolio, R_p, over [t - Δ t,t] is  

    R_p(t) = V^h(t) - V^h(t-Δ t)/V^h(t - Δ t) = ∑_i=1^2h_i(t)R_i(t),

where R_i is the rate of return of individual assets. The log return of our portfolio, r_p is given by

    r_p(t) = h_1r_1(t) + h_2r_2(t),

where r_i(t) ≈ R_i(t), as we showed in equation (<ref>).

Denote by V^h(t) the value of the portfolio corresponding to the admissible strategy h∈𝒜^1. Then:

  (i) The expectation of portfolio return over [t- Δ t, t] is 

    E(r_p(t)) = h_1E[r_X(t)] + h_2E[r_Y(t)].

  (ii) The variance of portfolio return over [t - Δ ,t] is

    Var(r_p(t))    =    h_1^2Var[r_X(t)] + h_2^2Var[r_Y(t)] + 2h_1h_2Cov[r_Xr_Y(t)],


where r(X_t)=ln( X_t/X_t-Δ t) and r(Y_t)=ln( Y_t/Y_t-Δ t) the daily log returns of assets X and Y and

  * E(r_X(t)) = (μ - σ^2/2)Δ t is the expected return of the asset price X  over the horizon [t - Δ t,t];
    
  * E(r_Y(t)) = [ ln( ae^μΔ t + (Y_0 -a)e^-κΔ t) -  ce^(2μ+σ^2)Δ t  + de^(μ -κ +σηρ)Δ t/2(ae^μΔ t + (Y_0 - a)e^-κΔ t)^2
- ln(Y_t-Δ t)] -
(Y_0^2-c-d)e^2(η^2 - κ)Δ t/2(ae^μ t + (Y_0 - a)e^-κΔ t)^2 + 1/2 is the expected return of the asset price Y  over the horizon [t - Δ t,t];
    
  * Var(r_X(t)) = σ^2Δ t is the variance of return of asset price X over the horizon [t -Δ t, t];
    
  * Var(r_Y(t)) =  ce^(2μ+σ^2)Δ t/(ae^μΔ t + (Y_0 - a)e^-κΔ t)^2 + de^(μ+σηρ -κ)Δ t/(ae^μΔ t + (Y_0 - a)e^-κΔ t)^2
 +(Y_0^2-c-d)e^2(η^2 - κ)Δ t/(ae^μΔ t + (Y_0 - a)e^-κΔ t)^2 - 1 is the variance of return of asset price Y over the horizon [t -Δ t, t];
    
  * Cov(r_X(t)r_Y(t)) = ln(be^(μ+σ^2)Δ t + (X_0Y_0 - b)e^(σηρ - κ)Δ t/aX_0e^2μΔ t + (X_0Y_0- aX_0)e^(μ -κ)Δ t) is the covariance of returns of two asset prices X and Y over the horizon [t -Δ t, t].

See Appendix <ref>
The optimal weights for mean-variance criterion were derived in <cit.>. We state the following proposition from <cit.> applied to the cointelation model (<ref>).

The optimal solution for the problem in (<ref>) for cointelation model (<ref>) is:

    h^*(t) = 1/e'Σ^-1(t)eΣ^-1(t)e + τ[ Σ^-1(t)M(t) - e'Σ^-1(t)M(t)/e'Σ^-1(t)eΣ^-1(t)e ],

with e' = [1,1,1], M(t)=[E(r_X(t)), E(r_Y(t))] and covariance matrix is

    Σ(t) = 
      [        Var(r_X(t)) Cov(r_X(t),r_Y(t)); Cov(r_X(t),r_Y(t))        Var(r_Y(t)) ],

where the expressions for E(r_X(t)), E(r_Y(t)) and  Var(r_X(t)), Var(r_Y(t)), Cov[r_X(t),r_Y(t)] are given above.

Replacing these formulas for expectation, variance and covariance of the returns of asset prices in equation (<ref>), we get optimal strategies for mean-variance optimization problem. We will present numerical examples in Section <ref>.



 §.§ Stochastic control for pairs trading 


  §.§.§ Power utility maximization problem

We now use a stochastic control approach to the power utility maximization problem. Here we mainly follow <cit.>, but with modified dynamics for asset prices. More specifically, they assume the price dynamics of one of the assets is a geometric Brownian motion and model the log-spread as an Ornstein-Uhlenbeck process. We, however, assume the dynamics of asset prices are governed by the cointelation model in equation (<ref>), where one of the assets follow the geometric Brownian motion and the second asset mean reverts around the first one. 

Let (Ω, ℱ, P) be a complete probability space with a filtration (ℱ_t)_t≥ 0 generated by two-dimensional Brownian motion: (W, W̃). We consider the same market as in Subsection <ref>: two assets which follow the cointelation model (<ref>).

We assume an initial wealth v_0 > 0 at time t=0. Initial wealth is held in a margin account. For simplicity we assume that the interest rate for margin account is 0, r=0.  Margin account restricts how much one can short or long. The holdings are allowed to be adjusted continuously up to a fixed horizon T. The investment behavior is modelled by an investment strategy π = (π_1, π_2). Here, π_i(t), i = 1, 2, denotes the percentage of total wealth invested in i-th asset at time t (see equation (<ref>)). Let π_1(t), π_2(t) be respectively the portfolio weights for assets X and Y at time t. We only allow pairs trading: short one of the asset and long the other in equal dollar amount, i.e. π_1(t) = -π_2(t). In addition, we restrict our considerations to self-financing strategies. 


We define admissible control and controlled process as in <cit.>.

Given a subset U of ℝ^2, we denote by 𝒰_0 the set of all progressively measurable processes π = {π_t,t≥ 0} valued in U. The elements of 𝒰_0 are called control processes.

Denote by V^π(t) the value of portfolio corresponding to strategy π at time t, which is given by 

    V^π(t) = π_1(t)V^π(t)/X(t)X(t) + π_2(t)V^π(t)/Y(t)Y(t).

The dynamics of the portfolio value V^π associated with strategy π = (π_1, π_2) is given by

    dV^π(t)=V^π(t)[ π_1(t)dX(t)/X(t) + π_2(t)dY(t)/Y(t)]

Replacing the dynamics for X(t) and Y(t) into (<ref>) we get:

    dV^π(t) = V^π(t)[ π_1(μ dt +σ dW(t)) - π_1(κ(X(t)/Y(t) -1)dt + η dW̃(t))].

Denote Z(t):=X(t)/Y(t). For the cointelation model (<ref>) we obtain that Z(t) has the dynamics

    dZ(t) = [μ +η^2 -σηρ -κ(Z(t) -1)]Z(t)dt + Z(t)(σ dW(t) +η dW̃(t)).

By Ito's quotient rule:

    d(X(t)/Y(t)) = dX(t)/X(t)X(t)/Y(t) - dY(t)/Y(t)X(t)/Y(t) +  d⟨ Y,Y⟩_t/Y(t)^2X(t)/Y(t) - d⟨ X,Y⟩_t/X(t)Y(t)X(t)/Y(t).

Writing this in terms of Z(t) gives

    dZ(t) = Z(t)(μ dt + σ  dW_t - κ (Z(t) - 1)dt - η dW̃(t) +  η^2Y^2(t)/Y^2(t)dt - σηρ dt)      
     
    = [μ + η^2 - σηρ - κ(Z(t) - 1)]Z(t)dt + (σ dW(t) -η dW̃(t))Z(t),

which proves the lemma.


For each control process π∈𝒰_0 we rewrite the dynamics of two-dimensional state process, P=(V^π,Z), as follows

    dP(t) = a(t,P(t),π(t))dt + b(t,P(t),π(t))dB(t).

with initial value of P(t_0) = p_0 and B = (W,W̃) being the two-dimensional Brownian motion. The process P is called the controlled process.
Let [t_0,T] with 0≤ t_0<T<∞ be the relevant time interval and define Q := [t_0,T)×ℝ^2. The coefficient functions 

    a      :     Q × U →ℝ^2, 
    
    b      :     Q × U →ℝ^2× 2,

are all continuous. Further, for all π∈ U let a(·, ·, π) and b(·, ·, π) be in C^1(Q). We then define

Denoting 𝒜^2 the set of all admissible controls, we say a control {π(t)}_t∈[t_0,T] will be called admissible if the following conditions hold

  (i)∀ k ∈ℕ the integrability condition

    E( ∫_t_0^T |π(s)|^kds ) < ∞

is satisfied,

  (ii) the corresponding state process P^π satisfies

    E^t_0,p_0( sup_t ∈ [t_0,T]|P^π(t)|^k) < ∞,

  (iii) only pairs trading is allowed: short one of the asset and long the other

    π_1 = -π_2.

Since we consider self-financing portfolio, then by equation (<ref>) the dynamics of the state process, P = (V^π, Z), becomes

    dV^π(t)=V^π(t)[ (π_1[μ - κ(Z(t) -1)])dt  + π_1[σ dW(t)) + η dW̃(t)] ],      V^π(0)=v_0,      
     
    dZ(t)=[μ +η^2 -σηρ -κ(Z(t) -1)]Z(t)dt + [σ dW(t) - η dW̃(t)]Z(t),   Z(0)=z_0.


  §.§.§ Optimal investment strategy

We assume that an investor's preference is represented by the power utility function

    U(x) = 1/γ x^γ ,

with x ≥ 0 and risk aversion parameter γ < 1. Our aim is to maximize the objective functional J over all admissible controls, i.e. determine an admissible control π(·) such that for each initial value (t_0, v_0) the utility functional below is maximized:

    J(t_0,v_0, z_0;π) := 𝔼[U(V^π(T)) | V_t_0 = v_0, Z_t_0=z_0 ].

The optimization problem is to find ṽ(t,v,z) and π∈𝒜^2 such that

    ṽ(t,v,z) := sup_π(·) ∈𝒜^2 J(t,v,z,π) = J(t,v,z,π^*).

Consider the function G(t,v,z) such that G ∈ C^1,2(Q). The Hamilton-Jacobi-Bellman (HJB) equation corresponding to the stochastic control problem (<ref>) is 

    ∂ G/∂ t(t,v,z) + sup_π∈𝒜^2ℒ^πG(t,v,z) = 0,

subject to terminal condition 

    G(T,v,z) = v^γ.

The infinitesimal generator, ℒ^π G(t,v,z) in (<ref>) associated with the two dimensional state process P = (V,Z) is given by

    ℒ^π G(t,v,z)    =   1/2[π^2_1(σ^2 - 2σηρ +η^2) v^2 G_vv + 2π_1 (σ^2 - 2σηρ +η^2) vzG_vz +  
           (σ^2 - 2σηρ  + η^2)z^2G_zz] + [π_1[μ - κ(z-1)]]vG_v + 
          [ μ +η^2 -σηρ - κ(z - 1)]zG_z.

If there exists an optimal control π^*(·) then G coincides with the value function:

    G(t,v,s) = ṽ(t,v,z) = J(t,v;π^*).

Using separation ansatz we reduce a 3-dimensional HJB equation in (<ref>) to the following 2-dimensional PDE:

    σ̃(γ-1) ff_t  - 1/2σ̃^2γ z^2f_z^2 - 1/2γ[μ -κ(z-1)]^2f + 1/2σ̃(γ-1)z^2ff_zz - 
          σ̃γ [μ -κ(z-1)]zff_z + σ̃(γ -1)[μ + η^2 -σηρ -κ(z-1)]ff_z,
                                  with   f(T,z) = 1,   (t,z)∈ [0,T]×ℝ,    ∀ z ∈ℝ,

where σ̃ =  σ^2 - 2σηρ +η^2. 

The issue at this stage is that this PDE does not have a closed for solution. This is a non standard PDE, which is not high dimensional but is nonlinear which makes using finite difference methods or any standard numerical methods inadequate. For this reason we propose to use the "Deep Galerkin Method" to  solve the PDE in (<ref>). Once the solution is found, we can write the optimal strategy as

    π_1^*    =    - σ̃zG_vz + [μ - κ(z-1)]G_v/σ̃vG_vv 
     = - σ̃z(f_zv^γ -1γ) + [μ - κ(z-1)](fv^γ -1γ)/σ̃v(fv^γ -2γ (γ -1))
       =    - σ̃zf_z + [μ - κ(z-1)]f/σ̃f (γ -1) = - zf_z/(γ -1)f - [μ - κ(z-1)]/σ̃(γ-1).

See Appendix <ref> for the details.






 §.§ Deep learning for solving PDE in stochastic control 

Without an analytical solution to the non-standard 2-dimensional PDE in (<ref>), we approximate the solution with an algorithm "Deep Galekin Method" (DGM) proposed in <cit.>. DGM is a merger of the Galerkin method and deep neural network machine learning algorithm. The Galerkin method is a popular numerical method which seeks a reduced-form solution to a PDE as a linear combination of basis functions. The deep learning algorithm, or DGM, uses a deep neural network instead of a linear combination of basis functions. The algorithm is trained on batches of randomly sampled time and  space points, therefore it is mesh free.


  §.§.§ Brief review of DGM

In general case, consider a PDE with d spatial dimensions:

    ∂ u/∂ t(t,x;θ) + ℒu(t,x)= 0,    (t,x) ∈ [0,T]×Ω , 
           u(t,x) = g(t,x),    x ∈∂Ω , 
           u(t=0,x) = u_0(x),    x ∈Ω

where x ∈Ω⊂ℝ^d and ℒ is an operator of all the other partial derivatives. The goal is to approximate the U(t,x) with deep neural network f(t,x;θ). Here θ∈ℝ^K are the neural network parameters. We want to minimize the objective function associated to the problem (<ref>) which consists of three parts:

  1.  A measure of how well the approximation satisfies the PDE:

    ∂ f/∂ t(t,x;θ) - ℒf(t,x;θ)^2_[0,T]×Ω, ν_1.

  2. A measure of how well the approximation satisfies the boundary condition:

    ∂ f/∂ t(t,x;θ) - g(t,x)^2_[0,T]×∂Ω, ν_2.

  3. A measure of how well the approximation satisfies the initial condition:

    ∂ f/∂ t(0,x;θ) - u(0,x)^2_Ω, ν_3.

Here all three errors are measured in terms of L^2-norm, i.e. f(y)^2_𝒴,ν = ∫_𝒴|f(y)|^2ν(y)dy with ν(y) being a density on region 𝒴. 

The sum of all three terms above gives us the objective function associated with the training of the neural network: 

    J(f)    =   ∂ f/∂ t(t,x;θ) - ℒf(t,x;θ)^2_[0,T]×Ω, ν_1 + 
          ∂ f/∂ t(t,x;θ) - g(t,x)^2_[0,T]×∂Ω, ν_2 + ∂ f/∂ t(0,x;θ) - u(0,x)^2_Ω, ν_3.

Thus, the goal is to find a set of parameters θ such that the function f(t,x;θ) minimizes the error J(f). When the dimension d is large, estimating θ by directly minimizing J(f) is infeasible. Therefore, one can minimize the error J(f) using a machine learning approach: stochastic gradient descent, where we use a sequence of time and space points drawn randomly.
The algorithm for DGM method is described in Algorithm <ref> below.
indent=2em
The learning rate, α_n, is a configurable hyperparameter[In machine learning, a hyperparameter is a parameter whose value is set before the learning process begins whereas, the values of other parameters are derived via training.] used in the training of neural networks that controls how much to change the model in response to the estimated error. Each time the model weights are updated. Learning rate has a small positive value, often in the range between 0.0 and 1.0. Similar to <cit.>, we set α_0 = 0.001. Note that our learning rate α_n must decrease with n, see <cit.>, and a simple enough way to do that is by using an exponential weighted method where α_n α_n-1*λ with λ∈]0,1[.

The neural network (NN) architecture used in DGM is like a long short-term networks (LSTMs) though with small differences, see <cit.>. We describe below the architecture of this NN:

    S^1 = σ(𝐰^1·𝐱+b^1) 
       Z^l = σ(𝐮^z,l·𝐱+𝐰^z,l· S^l+b^z,l)     l = 1, …, L
       G^l = σ(𝐮^g,l·𝐱+𝐰^g,l· S^l+b^g,l)     l = 1, …, L 
       R^l = σ(𝐮^r,l·𝐱+𝐰^r,l· S^l+b^r,l)     l = 1, …, L
       H^l = σ(𝐮^h,l·𝐱+𝐰^h,l· (S^l⊙ R^l)+b^h,l)     l = 1, …, L 
       S^l+1 = (1-G^l)⊙ H^l + Z^l ⊙ S^l     l = 1, …, L
       f(t,𝐱,θ) = 𝐰· S^L+1 + 𝐛

with ⊙ denoting Hadamard multiplication, L number of layers and σ the activation function. The rest of the subscript refer to the neurones for our NN architecture of Figures <ref> and <ref>.

We can see the Bird Eye view of the DGM <cit.> method in Figure <ref> and its details in Figure <ref>. The rationale is explained in <cit.>.


  §.§.§ Testing DGM on Merton problem

The method was tested with several nonlinear, high-dimensional PDEs independently in <cit.> and <cit.>, including nonlinear HJB equations. We have tested the DGM algorithm on HJB equation for the Merton problem ourselves. More specifically, Figures <ref> and <ref> show the plots of the analytical and approximated surface with DGM solution. Figure <ref> shows the difference between analytical and approximate solution. The approximation is good. Most of the time, the error is between 0% and 1%. The approximate solution does not do as well around t=0 (the maximum error of 4% is around t=0). This corroborates with the findings in <cit.>.



















  §.§.§ Solution to our PDE problem using DGM

Recall the PDE we want to solve is given in equation (<ref>). In the absence of a closed form solution to this PDE we approximate the solution with the DGM algorithm described above. Figure <ref> shows the approximate solution to the PDE in (<ref>) for different parameter values. Recall, once we have the numerical solution f for the PDE above, we obtain the optimal weights as following:

    π_1^*    =    - σ̃zG_vz + [μ - κ(z-1)]G_v/σ̃vG_vv = - σ̃z(f_zv^γ -1γ) - [μ - κ(z-1)](fv^γ -1γ)/σ̃v(fv^γ -2γ (γ -1))
       =    - σ̃zf_z - [μ - κ(z-1)]f/σ̃f (γ -1) = - zf_z/(γ -1)f - [μ - κ(z-1)]/σ̃(γ-1),

with π_1^* = - π_2^*.



 §.§ Dynamic Switching between optimal strategies of mean-variance and power utility

Although in the previous two cases we assume that an investor has a certain risk preferences as modelled by a utility function (MVC and power utility), it is interesting to consider a limiting case where the investor can be always persuaded to go for more money (identical utility function U(x) = x, which is essentially the power utility function with risk aversion parameter γ =1) when deciding between MVC or power utility.

Assuming that an investors' preference  is modelled either as in equation (<ref>) or as in equation (<ref>), in order to improve further the portfolio returns we employ dynamic switching between the two optimal strategies

    ψ^*(t) = π^*(t),   if  V^π^*(t) ≥ V^h^*(t), 
    
    h^*(t),   otherwise,



where π^*(t) and h^*(t) are given in equations (<ref>) and (<ref>) and V^π^*(t) and V^h^*(t) are given in equations (<ref>) and (<ref>). The motivation behind the dynamic switching is that the investor wants to benefit from both the mean-reversion and the correlation elements of the cointelation model (<ref>). More specifically, as the spread between two assets increases the investor implements pairs trading and makes profit, otherwise the MVC approach is used. 

The portfolio return over investment horizon [0,T] with T=1000 days is

    R(r_p) = V(0) - V(T)/V(0).

We perform 500 simulations with the same model and present in Table <ref> the average results. The average return at terminal time T obtained by using dynamic switching optimal strategies is higher than the average returns calculated by employing MVC or power utility maximizing optimal strategies.



§ MACHINE LEARNING FORMULATION OF THE PORTFOLIO OPTIMIZATION PROBLEM


 §.§ The portfolio optimization problem

We assume an initial wealth w̃_0 > 0 at time t=0. The investment behaviour is modeled by an investment strategy w = (w_1, w_2). Here w_1(t), w_2(t) denote the percentages of wealth invested in asset X and Y respectively at time t. Let V(t) denote the portfolio value at time t and V^PnL(t) := V(t) - V(0) denote the profit ant loss (P&L) over [0,t]. At each time t we allow either pairs trading: w_1(t) = -w_2(t) or long only strategies without leverage: w_1(t) + w_2(t) = 1 with w_1(t), w_2(t) > 0.

The general optimization problem is to find an optimal strategy, w(t), such that the
terminal P&L is maximized:  

    w^*(t) = w(t) ∈𝒜argmax V^PnL(w,T),

where V^PnL(w,T) is profit and loss corresponding to the strategy w at time terminal time T. We use clustering analysis to device the bands and in each band we solve the following optimization problem

    w_i^*(t) = w_i(t) ∈𝒜argmax V^PnL(w_i,t),

where i=1,...,n is the number of bands, V^PnL(w_i,t) is  profit and loss corresponding to the strategy w_i at time t. Then the overall solution w^* is obtained via a linear interpolation of optimal weights per band w_i^*

The advantage of the proposed method is that we do not impose certain model on the asset prices. Only data observations are required to calculate the optimal weight, meaning that the complex SDE calibration is avoided.



 §.§ Review of Band-Wise Gaussian Mixture model 

We review band-wise Gaussian mixture model because it inspires our method of selecting the bands. Consider a probability space (Ω, ℱ, ℙ) and let (P_t)_t ≥ 0 denote the asset price. Mahdavi-Damghani and Roberts <cit.> has recently introduced a generalised bumping SDE for the price dynamics of asset P_t. The SDE contains some secondary parameters whose purpose is empirical manual fitting. The generalized SDE is given by

    dP_t= θ_t,τ (μ_t,τ - P_t) dt + σ P_t^α(1-P_t^2)^β dW_t.

Here θ_t is the speed of mean reversion, μ_t is the long term mean, α is the positivity flag enforcer, β is the [-1,+1] boundary flag enforcer and {⋃ dW_i}_i=t-τ^t is the set of historical deviations of the assumed model's distribution (e.g.: all the historical absolute returns in the context of a normal diffusion).

This generalised SDE gives as a special case the cointelation model: take θ = -μ, μ = 0, α = 1 and β = 0 for the dynamic of X in (<ref>); take θ_t, τ = κ, μ_t,τ = X_t, α = 1, β = 0 for the dynamics of Y in (<ref>). The SDE in (<ref>) can also model:

  *  Proportional returns (log-normal diffusion) when θ=0, α=1, β=0.
	
  *  Absolute returns (normal diffusion) when θ=0, α=0, β=0,
	
  *  Mean reverting returns where we enforce positivity of returns (e.g. CIR <cit.> diffusion when α = 1/2 and β = 0),
	
  *  Mean reverting returns where we do not enforce positivity of the returns (e.g OU <cit.> diffusion when α = 0 and β = 0).

In general calibrating parameters of the SDE in  (<ref>) to a real data is complex. Using data simulated with (<ref>) their empirical distribution is approximated for the purpose of prediction by a band-wise Gaussian mixture model. This is done for a sequence of bands which are created using Machine Learning clustering method (see <cit.>).

Let P={p_1, …, p_n} be a set of empirical random variables sampled using equation (<ref>) with cumulative distribution function F(p) and density f(p). Denote O={p_(1), …, p_(n)} the ordered set of P such that p_(1) < p_(2) < …  < p_(n) and 

    O^i_h={p_(n((i-1)+1)/h), …, p_(n(i)/h)}.

Then the band-wise Gaussian mixture model for the empirical distribution function of the data simulated using the SDE in equation (<ref>) is given as follows:

    F̂_n(p_i|ℱ_t)=1/n∑ _j=1^h∑ _i=η^ζ1_p_i∈ O^j_h

with η = n((i-1)+1)/h and ζ=n(i)/h.

For example in the case bands h=3, using a Gaussian Mixture such that

    F̂_n(p_i|ℱ_t)=𝒩(-3,1)1_p_t∈ O^1_3+𝒩(0,1)1_p_t∈ O^2_3+𝒩(3,1)1_p_t∈ O^3_3,

we obtain the approximate stratification in Figure <ref>. The stratification is made so that the cardinality in each O^j_h region remains approximately the same, as opposed to being the result of a geometrical separation function of p_(1) and p_(n).

Theorem 1 in <cit.> ensures a good approximation of the generalised SDE (<ref>) by the Gaussian mixture model (<ref>). The calibration for the band-wise Gaussian mixture is given in Algorithm 
<ref>.

For our optimization problem we take a similar approach of dividing the range of observations into bands via the clustering algorithm, and then perform an optimization in each band via perturbation of weights.
indent=2em

 §.§ Optimal Machine Learning strategy

Based on the idea of band-wise Gaussian mixture model, we use clustering analysis to create bands, however not for the observed asset price data, but for the spread between two asset prices in (<ref>), i.e. X_t - Y_t. Inside of each band instead of specifying the distribution as in band-wise Gaussian mixture, we test a set of strategies that maximizes the corresponding P&L. We record the optimal strategies within each band, and in live trading, whenever the spread of asset prices falls in a certain band we employ the optimal strategy for this specific band.

We now present the trading signal that translates to investment strategy in machine learning approach.


  
*The Bayesian set-up:
We set from equation (<ref>)B_t = X_t - Y_t and have

    B_t = {B^+_n,t, B^+_n-1,t, …, B^+_1,t, B^1_1,t, …, B^-_n-1,t, B^-_n,t},

such that B^+_n,t > B^+_n-1,t > … > B^+_1,t > 0 > B^-_1,t > … > B^-_n-1,t > B^-_n,t.
We know that depending on the spread, the resulting approximated distribution of the samples differ <cit.>. The calibration algorithm will then consist of creating as many zones as possible whilst and as many strategies as possible within these bands and test how well each strategy is doing in each band in terms of P&L maximization. We take a direct approach (see Remark <ref>) consisting of 3 strategies and their cumulative P&Ls. Fixing the bands [a_i, b_i], with i =1, 2, ...,n we consider the following strategies:

  *  Strategy S^++ in which we are long both X and Y at time t in between bands [a_i,b_i] and with P&L V^++_[a_i,b_i],t.

  *  Strategy S^+- in which we are long X and short Y at time t in between bands [a_i,b_i] and with P&L V^+-_[a_i,b_i],t.

  *   Strategy S^-+ in which we are short X and long Y at time t in between bands [a_i,b_i] and with P&L V^-+_[a_i,b_i],t.

The P&Ls corresponding to these strategies are defined as following:

    V^++_[a_i,b_i],T   =   ∑_t=0^T [w^++_[a_i,b_i],tΔ X_t + (1-w^++_[a_i,b_i],t) Δ Y_t] 1_a_i< Δ_t≤ b_i, 
    
    V^+-_[a_i,b_i],T   =   ∑_t=0^T [w^+-_[a_i,b_i],tΔ X_t - (1-w^+-_[a_i,b_i],t) Δ Y_t] 1_a_i< Δ_t≤ b_i, 
    
    V^-+_[a_i,b_i],T   =   ∑_t=0^T [-w^-+_[a_i,b_i],tΔ X_t + (1-w^-+_[a_i,b_i],t) Δ Y_t] 1_a_i< Δ_t≤ b_i.

We call this approach direct, since ideally the number of strategies should consist of a more granular weight distribution. However for the sake of comparing with Financial Mathematics approach we consider the same set of strategies: long only, long/short.

We denote the maximum P&L achieved by each of these strategies by V^∓∓,*_[a_i,b_i],T, as given by equation (<ref>) and define S^**_[a_i,b_i],T of P&L V^**_[a_i,b_i],T (equation (<ref>)), the optimal strategy using Gaussian Learning in band [a_i,b_i]. 

    V^∓∓,*_[a_i,b_i],T = w^∓∓_[a_i,b_i],t∈[0,T]argmax V^∓_[a_i,b_i],T,  w^∓∓_[a_i,b_i],t∈ [0,1]

    V^**_[a_i,b_i],T = max (V^++,*_[a_i,b_i],T, V^+-,*_[a_i,b_i],T, V^-+,*_[a_i,b_i],T).

In live trading we recombine the optimal weights per bands into an overall optimal solution via a linear interpolation:

    w^*(t) = ∑_i=1^n w^*_i 1_{(X_t - Y_t) ∈ [a_i,b_i]}.

Although we do not have a proof that the resulting interpolated strategy in (<ref>) is optimal, we use it as a benchmark that still improves over the results with Financial Mathematics approach. Our goal is to apply Machine Learning approach to a pair of assets that exhibit some dependence, but this approach can be used for any model, i.e. it is model agnostic.
indent=2em

We further provide Algorithm <ref> as the pseudo-code for the calibration process. Note that in both Algorithms <ref> and <ref>, we have used a QuickSort which can be substituted by other sorting algorithms. Note that the use of self explanatory functions such as returnCorrespondingStrat(x,y) in line 20 of Algorithm <ref> which given the set of strategies and the P&L returns, as its name indicates, outputs the corresponding strategy that maximizes P&L. The function forecast(x,y,z) in line 22 of Algorithm <ref> takes as input the set of trained strategies and the current level of X_t and Y_t and returns a prediction of where the signals for the latter two should be. 

Finally the use of the argmax function in lines 13-16 can be replaced by a simple for loop but in the interest of not making the pseudocode too crowded we have kept it this way.


In <cit.> authors show that a reasonable risk manager or trader can assume the generalized SDE (<ref>) with β=0 and an α=1, in order to enforce positivity for the simulated scenarios of our risk factor. This very reasonable assumption would have crashed the whole risk engine if it is no longer satisfied in the real markets. The approach we advocate would have, however, been able to continue its dynamical learning scenario without any problem since it is model agnostic. 


§ NUMERICAL RESULTS 

Figure <ref> illustrates the ML and the DS approaches on one single simulated path. Note that when implementing the ML approach with a horizon of 1000 days, we double this data for training, i.e. we use 2000 historical daily prices. We have performed two sets of 500 simulations and we have gathered their results in the following two examples. 

 We have simulated 500 paths of X and Y based on cointelation model (<ref>) with parameters μ=0.05, σ=0.17, η=0.16, κ=0.1, ρ = -0.6. Figure <ref> illustrates that the Machine Learning approach with long/short strategies (ML_LS), on average performs slightly better in terms of P&L than the Stochastic Control approach (SC). However, based on histogram none of the approaches perform  significantly better or significantly worse than the other at any time. 
We have simulated 500 paths of X and Y based on cointelation model (<ref>) with parameters μ=0.05, σ=0.17, η=0.16, κ=0.1, ρ = -0.6. Figure <ref> illustrates how the ML approach seems to perform slightly better in terms of P&L than the FM approach about 55% of the time, while being outperformed the other 45% of the time. However, based on histogram we have noticed that sometimes the ML approach is being outperformed significantly more than it outperforms FM approach. 
 

From histogram of performance in Figures <ref> and <ref> we have concluded that for parameters  μ=0.05, σ=0.17, η=0.16, κ=0.1, ρ = -0.6 of cointelation model (<ref>) we have the following rankings for the approaches:

    SC < ML_LS < FM< ML.

The reason for ML with full set of strategies (long only and long/short) outperforming the DS most of the time might be the fact that in long only optimal strategies of ML approach we have more variety in weights, whereas the closed form formula (<ref>) in FM gives us almost constant weights (with small fluctuations).


 §.§ Possible directions for future work


  §.§.§ Multidimensional case
 
One direction for future work is to consider portfolio optimization problem for n-dimensional cointelation model. For instance, when n=3 we can have something of the following form:

    dS_t^a    =   σ S_t^a dW_t^a 
    
        dS_t^b    =   θ (S_t^a - S_t^b) dt + σ S_t^a dW_t^b 
     
        dS_t^c    =   θ (S_t^a - S_t^c) dt + σ S_t^c dW_t^c

One natural question would first be about how to model this triplet? For instance would equation (<ref>) with S^b and S^c reverting around S^a be more in-line with the pair from equation (<ref>) or would S^b reverting around S^a and S^c reverting around S^b be better? Are they equivalent or is one more useful? What happens as n increases? We plan to examine these questions in the future. 


























































  §.§.§ Application to cryptocurrencies

Another direction for future work is to model cryptocurrencies prices with the cointelation model and construct cryptocurrency indices using the portfolio optimization approaches proposed in this paper. Cryptocurrencies offer a source of alternative alpha, therefore there has been an emergence of cryptocurrency indices in the recent past with construction methodology ranging in the spectrum of Risk Parity to Stochastic Portfolio Theory (SPT). Given the spectacular volatility of the cryptocurrency market, even though the point of the index is to reduce the overall volatility, the index position remains fundamentally long. However, using a combination of beta neutral approach (long/short strategies) with an occasional long only alternative could be the winning combination for this asset class. For this reason the cointelation model is an interesting model to use for this application.



 §.§ Conclusion

We have studied the portfolio optimization problem of two assets that follow the cointelation model using two approaches: Financial Mathematics and Machine Learning. We first implemented the FM approach, where we use classic financial mathematics criteria: mean-variance and power utility maximization. Without an analytical solution to the PDE (<ref>), we resort to the DGM method, a deep learning algorithm, to solve it numerically.
The second approach we implemented is ML using clustering. The latter approach is easier to implement, it is model agnostic, therefore avoids the complex SDE calibration. In our case the Machine Learning approach slightly outperforms the Financial Mathematics approach.


§ APPENDICES


§ PROOF OF LEMMA <REF>

Since X_t is a geometric Brownian motion, we have 

    E[r(X_t)] = (μ - σ^2/2)Δ t

where X_t-Δ t is a known constant at time t - Δ t.
The expectation of log return of asset Y is

    E[r(Y_t)] = E[ln(Y_t)] - ln(Y_t-Δ t),

where Y_t-Δ t is a known constant at time t - Δ t. We use Taylor expansion to approximate expected value and variance of ln(Y_t) and covariance of ln(Y_t) and ln(X_t) (see <cit.>, p.165-167):

    E[ln(Y_t)]    ≈   ln(E[Y_t]) - σ^2[Y_t]/2E[Y_t]^2, 
    σ^2[ln(Y_t)]    ≈   σ^2[Y_t]/E[Y_t]^2,
    σ[ln(Y_t)ln(X_t)]    ≈   ln(1+σ[X_tY_t]/E[X_t]E[Y_t]).

First, we need to derive E[Y_t]. From equation (<ref>) we have

    Y_t = Y_t-Δ t + κ∫_t-Δ t^t(X_s - Y_s)ds + η∫_t-Δ t^tY_sdZ_s.

Taking expectation on both sides we have

    E[Y_t] = Y_t-Δ t + κ∫_t-Δ t^tE[X_s - Y_s]ds.

Differentiating on both sides we get

    dE[Y_t]/dt = κ E[X_t] - κ E[Y_t] = κ X_t-Δ t e^μΔ t -κ E[Y_t].

Denoting E[Y_t]=y(t) we obtain an ordinary differential equation (ODE):

    y' = -κ y + κ X_t-Δ te^μΔ t.

The solution is given by 

    y(t) = E[Y_t] = ae^μΔ t + (Y_t-Δ t - a)e^-κΔ t,

where 

    a = κ X_t-Δ t/μ+κ.

In order to derive E[Y_t^2] we first compute E[X_tY_t]. Applying integration by parts (IBP) to (<ref>) we get

    d(X_tY_t) = X_tdY_t + Y_tdX_t + dX_tdY_t 
           = κ X_t^2dt - κ X_tY_tdt + η X_tY_tdW_t +  μ X_tY_tdt +σ X_tY_tdZ_t + σηρ X_tY_tdt.

Thus 

    X_tY_t = X_t-Δ tY_t-Δ t + κ∫_t-Δ t^tX_s^2ds + η∫_t-Δ t^tX_sY_sdW_s + 
    
    (μ - κ +σηρ)∫_t-Δ t^t X_sY_sds + σ∫_t-Δ t^t X_sY_sdZ_s.

Taking expectation and differentiating on both sides 

    dE[X_tY_t]/dt = κ E[X_t^2] +(μ - κ + σηρ)E[X_tY_t].

Denoting E[X_tY_t]=x(t) we obtain ODE

    x' = κ E[X_t^2] + (μ + σηρ - κ)y.

Since X_t is GBM, its second moment is given by

    E[X_t^2] = E[X_t-Δ t^2e^(2μ - σ^2)Δ t + 2σ W_t] = X_t-Δ t^2e^(2μ+σ^2)Δ t.

Thus (<ref>) becomes

    x' = κ X_t-Δ t^2e^(2μ + σ^2)Δ t + (μ - κ + σηρ )y.

Using variation of parameters method we get the solution

    x(t) = E[X_tY_t] = be^(2μ+σ^2)Δ t + (X_t-Δ tY_t-Δ t - b)e^(μ - κ + σηρ )Δ t,

where 

    b = κ X_t-Δ t^2/μ+σ^2 + κ - σηρ.

Now we are ready to compute E[Y_t^2]. By Itô's lemma the dynamics of Y_t^2 is

    dY_t^2 = 2Y_tdY_t + (dY_t)^2 = (η^2 - 2κ)Y_t^2dt + 2κ X_tY_tdt + 2η Y_t^2dZ_t.

Integrating on both sides

    Y_t^2 = Y_0^2 + 2(η^2 -κ)∫_0^tY_s^2ds + 2κ∫_0^t X_sY_sds +2η∫_0^tY_s^2dZ_s.

Taking expectation on both sides and and differentiating 

    dE[Y_t^2]/dt =  2(η^2 - κ)E[Y_t^2] + 2κ E[X_tY_t].

Defining E[Y_t^2]=z(t) and replacing the value for E[X_tY_t] form equation (<ref>) we obtain an ODE

    z' = (η^2 - κ)z + 2κ be^(2μ+σ^2)Δ t + 2κ(X_t-Δ tY_t-Δ t - b)e^(μ - κ + σηρ )Δ t.

Using again variation of parameters we obtain the following solution

    z(t) = E[Y_t^2] = ce^(2μ+σ^2)Δ t + de^(μ -κ + σηρ )Δ t + (Y_t-Δ t^2-c-d)e^2(η^2 - κ)Δ t,

with c=2κ b /2μ+σ^2-2η^2 +2κ and d = 2κ(X_t-Δ tY_t-Δ t -b)/μ - 2η^2+κ+σηρ.

Now we are ready to approximate E[ln(Y_t)]. From (<ref>) we have

    E[ln(Y_t)] ≈ln[E[Y_t]] -E[Y_t^2]/2E[Y_t]^2+ 1/2 = ln( ae^μΔ t + (Y_t-Δ t -a)e^-κΔ t) + 1/2
           -  ce^(2μ+σ^2)Δ t + de^(μ -κ +σηρ)Δ t/2(ae^μΔ t + (Y_t-Δ t - a)e^-κΔ t)^2 -  (Y_t-Δ t^2-c-d)e^2(η^2 - κ)Δ t/2(ae^μΔ t + (Y_t-Δ t - a)e^-κΔ t)^2

and 

    E[r(Y_t)] = E[ln(Y_t)] - ln(Y_t-Δ t) ≈ln( ae^μΔ t + (Y_t-Δ t -a)e^-κΔ t) + 1/2 -
           ce^(2μ+σ^2)Δ t + de^(μ -κ +σηρ)Δ t/2(ae^μΔ t + (Y_t-Δ t - a)e^-κΔ t)^2 - (Y_t-Δ t^2-c-d)e^2(η^2 - κ)Δ t/2(ae^μΔ t + (Y_t-Δ t - a)e^-κΔ t)^2 - ln(Y_t-Δ t)

From (<ref>) we have
 
    Var[r(Y_t)] = Var[ln(Y_t)] ≈E[Y_t^2]/E[Y_t]^2 - 1 = 
           ce^(2μ+σ^2)Δ t + de^(μ+σηρ -κ)Δ t/(ae^μΔ t + (Y_t-Δ t - a)e^-κΔ t)^2 + (Y_t-Δ t^2-c-d)e^2(η^2 - κ)Δ t/(ae^μΔ t + (Y_t-Δ t - a)e^-κΔ t)^2 -1

and  
 
    Var[r(X_t)] = Var[ln(X_t)] = σ^2 Δ t.

From (<ref>) we obtain the covariance:

    Cov[r(X_t)r(Y_t)] = Cov[ln(X_t)ln(Y_t)] ≈ln(E[X_tY_t]/E[X_t]E[Y_t])  
       ≈ln(be^(2μ+σ^2)Δ t + (X_t-Δ tY_t-Δ t - b)e^(μ -κ + σηρ)Δ t/aX_0e^2μΔ t + (Y_t-Δ tX_t-Δ t- aX_t-Δ t)e^(μ-κ)Δ t).


§ DIMENSION REDUCTION OF 3-DIM HJB (<REF>)
 For ease of notation let σ̃ = σ^2 - 2σηρ +η^2 and rewrite (<ref>):

    G_t + sup_π_1{1/2(π_1^2 σ̃ v^2G_vv + σ̃z^2G_zz + 2π_1σ̃vzG_vz)+ 
          (π_1[μ - κ(z-1)])vG_v + (μ + η^2  - σηρ - κ(z-1))zG_z} =0.

The first order condition for the maximization is

    π_1^*σ̃vG_vv + σ̃zG_vz + [μ - κ(z-1)]G_v = 0.

Now assuming G_vv<0 the first order condition is sufficient, yielding

    π_1^* = - σ̃zG_vz + [μ - κ(z-1)]G_v/σ̃vG_vv.

Replacing (<ref>) back into (<ref>) yields:

    G_t + 1/2{(σ̃zG_vz + [μ - κ(z-1)]G_v)^2/σ̃^2v^2G^2_vvσ̃v^2G_vv + σ̃z^2G_zz
          - 2σ̃zG_vz + [μ - κ(z-1)]G_v/σ̃vG_vvσ̃vzG_vz} + [μ +η^2 -σηρ -κ(z-1)]zG_z 
          - σ̃zG_vz + [μ - κ(z-1)]G_v/σ̃vG_vv[μ - κ(z-1)]vG_v   = 0.

Multiplying both sides of equation by σ̃G_vv we get:

    σ̃G_tG_vv + 1/2(σ̃zG_vz + [μ -κ(z-1)]G_v)^2 - (σ̃zG_vz + [μ - κ(z-1)]G_v)σ̃zG_vz
          + 1/2σ̃z^2G_zzG_vv - (σ̃zG_vz + [μ - κ(z-1)]G_v)[μ -κ(z-1)]G_v 
          + σ̃[μ +η^2 -σηρ -κ(z-1)]zG_zG_vz = 0.

Expanding gives

    σ̃G_tG_vv + 1/2σ̃^2z^2G^2_vz + 1/2[μ -κ(z-1)]^2G^2_v + σ̃z[μ -κ(z-1)]G_vG_vz
            -σ̃^2z^2G^2_vz - σ̃z[μ -κ(z-1)]G_vG_vz +1/2σ̃z^2G_zzG_vv - σ̃z[μ -κ(z-1)]G_vG_vz
           - [μ -κ(z-1)]^2G^2_v + σ̃[μ +η^2 -σηρ -κ(z-1)]zG_zG_vz = 0.

Which further simplifies to

    σ̃G_tG_vv -1/2[μ -κ(z-1)]^2G^2_v + 1/2σ̃z^2G_zzG_vv - σ̃z[μ -κ(z-1)]G_vG_vz
           -  1/2σ̃^2z^2G^2_vz + σ̃[μ +η^2 -σηρ -κ(z-1)]zG_zG_vz = 0.

At this stage we were able to turn our four variable PDE into three, but we can get eliminate one more. For this we consider the following separation ansatz:

    G(t, v, z) = f(t,z)v^γ,

with the terminal condition

    f(T,z) = 1      ∀ z.

We compute the derivatives of (<ref>):

    G_t = f_tv^γ,    G_v = fv^γ -1γ ,   G_z=f_zv^γ ,   G_vv = fv^γ -2γ (γ -1), 
        G_vz = f_zv^γ -1γ ,    G_zz =v^γ f_zz.

Replace derivative back into (<ref>) and divide by v^2(γ -1)γ to get

    σ̃(γ-1) ff_t  - 1/2σ̃^2γ z^2f_z^2 - 1/2γ[μ -κ(z-1)]^2f +1/2σ̃(γ-1)z^2ff_zz -
          σ̃γ [μ -κ(z-1)]zff_z +  σ̃(γ -1)[μ + η^2 -σηρ -κ(z-1)]ff_z=0.
 
We now have a PDE with only two variables instead of four.





9SolvingPDEwithDL [1] Al-Aradi, A., Correia, A., Naiff, D., Jardim, G. and Sapito, Y. Solving nonlinear
and high-dimensional partial differential equations via deep learning. (2018) arXiv:1811.08782Alexander_Optim_heging [2] Alexander, C. Optimal hedging using cointegration. Philosophical Transactions of
the Royal Society a Mathematical Physical and Engineering Sciences 357 (1999), 2039–2058.

Alexander [3] Alexander, C.Market models: a guide to financial data analysis. John Wiley Sons, Inc., Chichester, West Sussex, (2001).

Benaroya [4] Benaroya, H., Han, S.M. and Nagurka, M. Probability models in engineering and science.  CRC Press, (2005).

Bodie [5] Bodie, Z., Kane, A. and Marcus, A. Investments. 4th edition, Irwin/McGraw-Hill, Chicago, (1999).

CIR [6] Cox, J.C., Ingersoll, J.E. and Ross, S.A. A theory of the term structure of interest rates. Econometrica 53 (1985), 385-407.

EngleGranger [7] Engle, R. and Granger, C. Long-run Economic Relationships: Readings in Cointegration. Oxford University Press, Oxford, New York, (1991).

Garciaa [8] Garcia, R., Gonzalez, V., Contreras, J. and Custodio, J. Applying modern portfolio theory for a dynamic energy portfolio allocation in electricity markets. Electric Power
Systems Research 150 (2017), 11–22.

HarrisonKreps [9] Harrison, J. and Kreps, D. Martingales and arbitrage in multiperiod securities markets. Journal of Economic Theory 20 (1979), 381–408.

KornKraft [10] Korn, R., and Kraft, H. A Stochastic Control Approach to Portfolio Problems with Stochastic Interest Rates. SIAM Journal on Control and Optimization 4 (2002), 1250-1269.

LiNg [11] Li, D. and Ng, W. Optimal dynamic portfolio selection: Multiperiod mean-variance formulation. Mathematical Finance 10  (2000), 387–406.

CointAndInferred [12] Mahdavi-Damghani, B. The non-misleading value of inferred correlation: An introduction to the cointelation model. Wilmott Magazine 67 (2013), 0-61.

bandwiseGaussianMix [13] Mahdavi-Damghani, B. and Roberts, S. A Proposed Risk Modeling Shift from the Approach of Stochastic Differential Equation towards Machine Learning Clustering: Illustration with the Concepts of Anticipative & Responsible VaR. (2017), SSRN:3039179.

Damgani [14] Mahdavi-Damghani, B.,  Welch, D., O'Malley, K. and Knights, S. The misleading value of measured correlation. Wilmott Magazine 61 (2013), 64-73.

MudchanatongsukPrimbsWong [15] Mudchanatongsuk, S., Primbs, J., and Wong, W. Optimal Pairs Trading: A Stochastic Control Approach. In American Control Conference (Seattle, Washington, USA, 2008), IEEE.

DGM [16] Sirignano, J. and Spiliopoulos, K. DGM: A deep learning algorithm for solving partial differential equations. Journal of Computational Physics 375 (2018), 1339 – 1364.

Soeryana [17] Soeryana, E., Fadhlina, N., Sukono, Rusyaman, E., and Supian, S. Mean-variance portfolio optimization by using time series approaches based on logarithmic utility function. In Materials Science and Engineering (2017), IOP.

OU [18] Uhlenbeck, G.E., and Ornstein, L.S. On the theory of Brownian Motion. Physical Review 35 (1930), 823-841.

Vidyamurthy [19] Vidyamurthy, G. Pairs Trading: Quantitative Methods and Analysis. John Wiley Sons, Inc., Hoboken, New Jersey, (2004).

Wilmott [20] Wilmott, P. Paul Wilmott Introduces Quantitative Finance. John Wiley Sons Ltd., Chichester, West Sussex, (2007).





  Figure <ref>: (Up) Simulated path of cointelation model (<ref>) with ρ=-1, θ=0.1, σ=0.01; (Down) Corresponding measured correlation (<ref>) as a function of the time increment increases from -1 to 1.
    
    
  Figure <ref>: Bird’s-eye perspective of overall DGM architecture <cit.>.
    
  Figure <ref>: Operations within a single DGM layer <cit.>
  Figure <ref>: Analytical solution of the Merton Problem.
    
  Figure <ref> Approximate solution of Merton problem using DGM.
    
  Figure <ref>: Error between analytical and approximate solution of Merton problem.
    
  Figure <ref>: Approximate solutions to PDE (<ref>) with DGM for four different scenarios of ρ and μ and fixed σ=0.2, η=0.19, γ=0.5.
    
  (a) Approximate solution with low μ = 0.01 and low ρ=-0.5.
        
  (b) Approximate solution with low μ = 0.01 and high ρ = 0.5.
        
  (c) Approximate solution with high μ = 0.4 and low ρ = -0.5.
        
  (d) Approximate solution with high μ = 0.4 and high ρ =0.5.
    
  Figure <ref>: Two examples of Gaussian Mixture Simulations with different number of bands. 
    
  (a) Empirical distribution of random variable sampled from cointelation model (<ref>) in three different zones described in Figure <ref>.
        
  (b) Empirical distribution of random variable sampled from cointelation model (<ref>) in five different zones: two additional zones were added to the initial three zones in Figure <ref>.
    
  Figure <ref>: (a) one simulated scenario based on cointelation model (<ref>) with parameters: μ=0.05, σ=0.17, η=0.16, κ=0.1, ρ = -0.6 and scaled spread: κ(X_t - Y_t); (b)  portfolio return and optimal weight of asset X with Dynamic Switching approach; (c) portfolio return and optimal weight of asset X and Y with Machine Learning approach.
    
  Figure <ref>: Histogram of excess (P&L) for ML_LS vs SC at terminal time T.
    
  Figure <ref>: Histogram of excess (P&L) for ML vs FM at terminal time T.
    
  Figure <ref>: In FM approach optimal long/short strategies are more volatile than optimal long only strategies.
    
  Figure <ref>: In ML approach optimal long/short strategies are slightly more volatile than optimal long only strategies.