#!/usr/bin/env python
# coding: utf-8

# # Portfolio Optimization
# 
# This notebook can be run online without installing any packages. Just click the logo:
# 
# [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/mcdeoliveira/pyoptimum-examples/master?filepath=examples%2Fportfolio.ipynb)
# 
# to run it on [binder](https://mybinder.org).
# 
# See this [case study](http://vicbee.net/portfolio.html) for more details.

# ## Import needed libraries

# In[1]:


import math
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd


# ## Initialize optimize.vicbee.net client

# In[2]:


username = 'demo@optimize.vicbee.net'
password = 'optimize'


# In[3]:


import pyoptimum # for easy access to the optimmize.vicbee.net api
client = pyoptimum.Client(username=username, password=password)


# # Classic Markowitz Portfolio Optimization
# 
# In the classic Markowitz optimal portfolio optimization (Markowitz, 1952), one seeks to determine 
# * $x$: a vector of porfolio positions,
# 
# that minimize the total portfolio variance given:
# * $r$: the vector of expected returns,
# * $Q$: the covariance matrix of the returns,
# * $\mu$: the desired portfolio return.
# 
# The problem can be formulated mathematically as the Quadratic Program (QP):
# $$
# \begin{aligned}
#   \min_{x} \quad & x^T Q x \\
#   \text{s.t.} \quad & r^T x = \mu \quad \sum_i{x_i} = 1 
# \end{aligned}
# $$
# 
# In this notebook we will show how one can setup and solve such portfolio optimization problems using the [Optimize API](https://optmize.vicbee.net/api/ui).
# 
# Since we are using Python, we will take advantage of [numpy](https://numpy.org/) and [pandas](https://pandas.pydata.org/) to manipulate the data in this notebook. We will also use [matplotlib](https://matplotlib.org/) to plot the data. The use of such libraries or even Python is not a requirement for using the API, which can be accessed using <a href="https://en.wikipedia.org/wiki/Ajax_(programming)">AJAX calls</a>. 
# 
# Data for the problem considered in this demo is taken from Elton and Gruber, 1995.
# 
# ### References:
# 
# 1. Harry M. Markowitz, “Portfolio Selection”. Journal of Finance, 7(1): 77 - 91, 1952.
# 2. Edwin J. Elton and Martin J. Gruber, Modern Portfolio Theory and Investment Analysis. Fifth Edition, John Wiley and Sons, 1995.
# 
# ## Setup data using pandas
# 
# Consider the following set of assets with associated expected return and return stadard deviation given as percentages assembled as a pandas DataFrame:

# In[4]:


data = [
    ['S & P',14.0,18.5],
    ['Bonds',6.5,5.0],
    ['Canadian',11.0,16.0],
    ['Japan',14.0,23.0],
    ['Emerging Markets',16.0,30.0],
    ['Pacific',18.0,26.0],
    ['Europe',12.0,20.0],
    ['Small Stocks',17.0,24.0] 
]
assets = pd.DataFrame(data, columns=['Label', 'Return (%)', 'STD (%)'])
assets


# The assets are associated with the following correlation coefficients:

# In[5]:


n = len(assets)
data = np.array([
    [1,.45,.7,.2,.64,.3,.61,.79],
    [.45,1,.27,-.01,.41,.01,.13,.28],
    [.7,.27,1,.14,.51,.29,.48,.59],
    [.2,-.01,.14,1,.25,.73,.56,.13],
    [.64,.41,.51,.25,1,.28,.61,.75],
    [.3,.01,.29,.73,.28,1,.54,.16],
    [.61,.13,.48,.56,.61,.54,1,.44],
    [.79,.28,.59,.13,.75,.16,.44,1]
])
cov = pd.DataFrame(data=data, columns=assets.Label.values, index=assets.Label.values)
cov


# ## Solve problem using the Optmize API
# 
# Before solving the problem we need to assemble the required data, that is the vector of expected returns $r$, which can be taken directly from the data:

# In[6]:


# vector of expected returns
r = assets['Return (%)'].values


# the covariance matrix $Q$, which can be calculated from the given data as follows:

# In[7]:


# covariance matrix
sigmas = np.diag(assets['STD (%)'].values)
Q = sigmas @ cov.values @ sigmas


# In the context of the classic Markowitz portfolio optimization we should also specify a desired expected return

# In[8]:


# expected return
mu = 14


# and then select a portfolio corresponding to this return with the smallest possible variance.
# 
# The actual calculation of the optimal portfolio will be done on the cloud.
# 
# 
# ## Submit problem to the api
# 
# We assemble all the data in an object and make an api call to get the optimal portfolio:

# In[9]:


# prepare data for submitting request to api
data = {
    'Q': Q.tolist(),
    'r': r.tolist(),
    'mu': mu
}
sol = client.call('portfolio', data)


# Once an optimal solution is calculated we append it to the dataframe.

# In[10]:


if sol['status']=='infeasible':
    print('> Problem is infeasible')
    assets['Optimal (%)'] = np.nan
else:
    print('> Problem is feasible')
    print('  * mu = {:.3f}%, std = {:.3f}% ({})'.format(mu, math.sqrt(sol['obj']), sol['status']))
    pd.options.display.float_format = '{:,.1f}'.format
    assets['Optimal (%)'] = 100*np.array(sol['x'])
assets


# By default, the solution vector adds to one, each number being interpreted as the fraction of the available budget that should be invested in order to realize the optimal portfolio. We use numpy to convert that to a percentage before displaying the solution.

# ## Solution without short sales

# The minimum variance portfolio that meets a required level of expected return is often one in which the positions can be negative. Negative positions are realized as *short sales*.
# 
# For example, in the problem solved above

# In[11]:


x = assets['Optimal (%)'];
print('{:.1f}%'.format(-x[x < 0].sum()))


# of the portfolio is comprised of short sales.
# 
# One can add constraints to the optimization problem in order to prevent short sales.
# 
# This can be done by adding the following option:

# In[12]:


# options
options = {
    'short': False
}


# and calling the API once again

# In[13]:


# prepare data for submitting request to api
data = {
    'Q': Q.tolist(),
    'r': r.tolist(),
    'mu': mu,
    'options': options
}
sol = client.call('portfolio', data)


# The new solution

# In[14]:


if sol['status']=='infeasible':
    print('> Problem is infeasible')
    assets['No Shorts (%)'] = np.nan
else:
    print('> Problem is feasible')
    print('  * mu = {:.3f}%, std = {:.3f}% ({})'.format(mu, math.sqrt(sol['obj']), sol['status']))
    assets['No Shorts (%)'] = 100*np.array(sol['x'])
assets


# is one without short sales:

# In[15]:


x = assets['No Shorts (%)']
print('{:.1f}%'.format(-x[x < 0].sum()))


# in which each every position is an asset purchase:

# In[16]:


print('{:.1f}%'.format(x.sum()))


# Note that the new portofolio is *riskier* than the optimal one with short sales.

# ## Accounting for cashflow

# Now suppose that one has an initial portfolio that they would like to rebalance, or buy or sell more assets.
# 
# In this case, one starts with an initial positions, for example:

# In[17]:


assets['Initial (%)'] = [20, 30, 0, 20, 0, 10, 0, 20]
assets


# An optimal portfolio will not change depending on whether one holds assets or not, so the answers obtained so far should not be affected by the existing holdings.
# 
# However, imagine that one has cash to invest and would like to know what is the best solution for investing the additional money without selling any of their assets.
# 
# One can take that scenario into account by adding an initial portfolio:

# In[18]:


# initial portfolio
x0 = assets['Initial (%)'].values/100


# and the desired *cashflow*, say equivalent to 40% of the existing portfolio,

# In[19]:


# cashflow
cashflow = 40/100


# along with the other data:

# In[20]:


# expected return
mu = 14

# options
options = {
    'sell': False 
}

# prepare data for submitting request to api
data = {
    'Q': Q.tolist(),
    'r': r.tolist(),
    'mu': mu,
    'x0': x0.tolist(),
    'cashflow': cashflow,
    'options': options
}
sol = client.call('portfolio', data)


# The optimal solution now is in the form of additional *buys* that need to be realized within the specified cashflow to achieve the desired return:

# In[21]:


if sol['status']=='infeasible':
    print('> Problem is infeasible')
    assets['Buys (%)'] = np.nan
    assets['After Cashflow (%)'] = np.nan
else:
    print('> Problem is feasible')
    print('  * mu = {:.3f}%, std = {:.3f}% ({})'.format(mu, math.sqrt(sol['obj']), sol['status']))
    x = np.array(sol['x'])
    assets['Buys (%)'] = 100*(x - x0)/cashflow
    assets['After Cashflow (%)'] = 100*x/x.sum()
assets


# ## Infeasible portfolio problems
# 
# As one add more constraints, it is often the case that certain problems become unsolvable.
# 
# For example, in the above cashflow problem, a solution does not exist if one lowers the ammount available for investment. For example, lowering the available cashflow to 10% or even 20% as in:

# In[22]:


# initial portfolio
x0 = assets['Initial (%)'].values/100

# cashflow
cashflow = 20/100

# expected return
mu = 14

# options
options = {
    'sell': False 
}

# prepare data for submitting request to api
data = {
    'Q': Q.tolist(),
    'r': r.tolist(),
    'mu': mu,
    'x0': x0.tolist(),
    'cashflow': cashflow,
    'options': options
}
sol = client.call('portfolio', data)


# results in an *infeasible* problem:

# In[23]:


if sol['status']=='infeasible':
    print('> Problem is infeasible')
    assets['Buys (%)'] = np.nan
    assets['After Cashflow (%)'] = np.nan
else:
    print('> Problem is feasible')
    print('  * mu = {:.3f}%, std = {:.3f}% ({})'.format(mu, math.sqrt(sol['obj']), sol['status']))
    x = np.array(sol['x'])
    assets['Buys (%)'] = 100*(x - x0)/cashflow
    assets['After Cashflow (%)'] = 100*x/x.sum()
assets


# Of course the reason why the problem has no solution is because the return of the initial portfolio:

# In[24]:


print('{:.1f}%'.format(r @ x0))


#  is much lower than the return asked of the optimization, which was of 14%.
#  
#  In this particular case, one should expect to obtain solutions with smaller cashflows if one asks for more modest returns.
#  
#  For example, a modest improvement in return from 12.75% to 13% is available with a 10% or 20% cashflow:

# In[25]:


# initial portfolio
x0 = assets['Initial (%)'].values/100

# cashflow
cashflow = 20/100

# lower return
mu = 13

# options
options = {
    'sell': False 
}

# prepare data for submitting request to api
data = {
    'Q': Q.tolist(),
    'r': r.tolist(),
    'mu': mu,
    'x0': x0.tolist(),
    'cashflow': cashflow,
    'options': options
}
sol = client.call('portfolio', data)


# The obtained optimal portfolio:

# In[26]:


if sol['status']=='infeasible':
    print('> Problem is infeasible')
    assets['Buys (%)'] = np.nan
    assets['After Cashflow (%)'] = np.nan
else:
    print('> Problem is feasible')
    print('  * mu = {:.3f}%, std = {:.3f}% ({})'.format(mu, math.sqrt(sol['obj']), sol['status']))
    x = np.array(sol['x'])
    assets['Buys (%)'] = 100*(x - x0)/cashflow
    assets['After Cashflow (%)'] = 100*x/x.sum()
assets


# is a step toward improving the returns with the available cashflow.
# 
# We will revisit this problem later in this notebook after showing how to calculate efficient frontiers.
# 
# ## General allocation constraints
# 
# It is also possible to specify arbitrary constraints on sets of assets.
# 
# For example, consider the that one desires an optimal portfolio with no cashflow, no shorts, but at least 50% of their assets on the combined S & P and Bond markets. Such a problem can be formulated after defining a set and corresponding constraint:

# In[27]:


constraint1 = {'label': 's&p_bonds', 'set': [0, 1], 'bounds': [.5, np.inf]}


# Note how sets are defined using numeric indices (starting from 0), and the constraint bounds estalish the lower and upper bound on the sum of the assets in the set.
# 
# The resulting constrained problem can be solved as before:

# In[28]:


# Expected return
mu = 14

# options
options = {
    'short': False 
}

# prepare data for submitting request to api
data = {
    'Q': Q.tolist(),
    'r': r.tolist(),
    'mu': mu,
    'constraints': [constraint1],
    'options': options
}
sol = client.call('portfolio', data)


# The obtained optimal portfolio:

# In[29]:


if sol['status']=='infeasible':
    print('> Problem is infeasible')
    assets['S&P + Bonds >= 50 (%)'] = np.nan
else:
    print('> Problem is feasible')
    print('  * mu = {:.3f}%, std = {:.3f}% ({})'.format(mu, math.sqrt(sol['obj']), sol['status']))
    x = np.array(sol['x'])
    assets['S&P + Bonds >= 50 (%)'] = 100*x
assets


# satisfies the constraint that

# In[30]:


print('{:.1f}%'.format(assets['S&P + Bonds >= 50 (%)'][0:2].sum()))


# is not less than 50%.
# 
# One might find it easier to formulate constraints based on the assets' labels. This is possible if a property label is added to the problem formulation:

# In[31]:


# Expected return
mu = 14

# options
options = {
    'short': False 
}

# prepare data for submitting request to api
data = {
    'Q': Q.tolist(),
    'r': r.tolist(),
    'mu': mu,
    'labels': assets['Label'].tolist(),
    'constraints': [{'label': 's&p_bonds', 'set': ['S & P', 'Bonds'], 'bounds': [.5, np.inf]}],
    'options': options
}
sol = client.call('portfolio', data)


# Note that no data is stored in the servers, and that labels can be easily anonymized if desired.
# 
# The above formulation produces the same optimal portfolio:

# In[32]:


if sol['status']=='infeasible':
    print('> Problem is infeasible')
    assets['S&P + Bonds >= 50 (%)'] = np.nan
else:
    print('> Problem is feasible')
    print('  * mu = {:.3f}%, std = {:.3f}% ({})'.format(mu, math.sqrt(sol['obj']), sol['status']))
    x = np.array(sol['x'])
    assets['S&P + Bonds >= 50 (%)'] = 100*x
assets


# satisfying the constraint

# In[33]:


print('{:.1f}%'.format(assets['S&P + Bonds >= 50 (%)'][0:2].sum()))


# ## Leverage constraints
# 
# Another type of constraint that can be easily imposed is leverage constraints.
# 
# We have already determined the optimal portfolio with and without short sales. The optimal portfolio has leverage at aroud 27%:

# In[34]:


print('{:.1f}%'.format(assets['Optimal (%)'][assets['Optimal (%)'] <= 0].sum()))


# while the optimal portfolio without short sales has, of course, 0% leverage.
# 
# How about an optimal intermediate portfolio with, say, 10% leverage? It can be determined by adding a *leverage constraint*:

# In[35]:


constraint = {'label': '10% leverage', 'function': 'leverage', 'set': ':', 'bounds': .1}


# Note the use of the function `leverage` and the set `:`, which is a shortcut representing all assets. The associated optimal porfolio is calculated below.

# In[36]:


# Expected return
mu = 14

# options
options = {
    'short': True 
}

# prepare data for submitting request to api
data = {
    'Q': Q.tolist(),
    'r': r.tolist(),
    'mu': mu,
    'labels': assets['Label'].tolist(),
    'constraints': [constraint],
    'options': options
}
sol = client.call('portfolio', data)


# In[37]:


if sol['status']=='infeasible':
    print('> Problem is infeasible')
    assets['Leverage <= 10 (%)'] = np.nan
else:
    print('> Problem is feasible')
    print('  * mu = {:.3f}%, std = {:.3f}% ({})'.format(mu, math.sqrt(sol['obj']), sol['status']))
    x = np.array(sol['x'])
    assets['Leverage <= 10 (%)'] = 100*x
assets


# which satisfies the desireded leverage requirements:

# In[38]:


print('{:.1f}%'.format(assets['Leverage <= 10 (%)'][assets['Leverage <= 10 (%)'] <= 0].sum()))


# ## Group leverage constraints
# 
# Leverage constraints can also be imposed in groups of assets.
# 
# The next example calculates an optimal portolio with at most 10% leverage in the Pacific and Europe regions.  

# In[39]:


# constraint
constraint = {'label': '10% group leverage', 'function': 'leverage', 'set': ['Pacific', 'Europe'], 'bounds': .1}

# Expected return
mu = 14

# options
options = {
    'short': True 
}

# prepare data for submitting request to api
data = {
    'Q': Q.tolist(),
    'r': r.tolist(),
    'mu': mu,
    'labels': assets['Label'].tolist(),
    'constraints': [constraint],
    'options': options
}
sol = client.call('portfolio', data)


# In[40]:


if sol['status']=='infeasible':
    print('> Problem is infeasible')
    assets['Group Levegare (%)'] = np.nan
else:
    print('> Problem is feasible')
    print('  * mu = {:.3f}%, std = {:.3f}% ({})'.format(mu, math.sqrt(sol['obj']), sol['status']))
    x = np.array(sol['x'])
    assets['Group Leverage (%)'] = 100*x
assets


# Note that leverage is calculated within the group's positions. In this case the group leverage is:

# In[41]:


group = assets['Group Leverage (%)'][5:7]
print('{:.1f}%'.format(100*group[group <= 0].sum()/group.sum()))


# as desired.
# 
# ## Sales constraints
# 
# Similarly to leverage constraints, it is possible to enforce constraints on the total amount of sales in a portfolio.
# Such problems are relevant in the context of cashflow problem in the presence of nonzero initial positions. With zero initial positions such sales constraints are equivalent to the leverage constraints introduced above.
# 
# For example, let us revisit the above cashflow problem for which we already know that a return of 14% is not possible with an injection of 20% of capital without selling assets. Next, we attempt to solve this same problem but this time imposing no short sales and a maximum of 25% in sales: 

# In[42]:


# constraint
constraint = {'label': '10% sales constraint', 'function': 'sales', 'set': ':', 'bounds': .25}

# initial portfolio
x0 = assets['Initial (%)'].values/100

# cashflow
cashflow = 20/100

# expected return
mu = 14

# options
options = {
    'sell': True,
    'short': False
}

# prepare data for submitting request to api
data = {
    'Q': Q.tolist(),
    'r': r.tolist(),
    'mu': mu,
    'x0': x0.tolist(),
    'cashflow': cashflow,
    'constraints': [constraint],
    'options': options
}
sol = client.call('portfolio', data)


# The optimal solution is:

# In[43]:


if sol['status']=='infeasible':
    print('> Problem is infeasible')
    assets['Trades (%)'] = np.nan
    assets['After Trades (%)'] = np.nan
else:
    print('> Problem is feasible')
    print('  * mu = {:.3f}%, std = {:.3f}% ({})'.format(mu, math.sqrt(sol['obj']), sol['status']))
    x = np.array(sol['x'])
    assets['Trades (%)'] = 100*(x - x0)/cashflow
    assets['After Trades (%)'] = 100*x/x.sum()
assets


# which satisfies both the no short sale constraints as well as the 25% cap on sales:

# In[44]:


trades = assets['Trades (%)']
print('{:.1f}%'.format(100*trades[trades <= 0].sum()/trades.sum()))


# # The efficient frontier
# 
# The efficient frontier is the set of minimal variance portfolios obtained after varying the expected returns. A graphical representation of the frontier can be ploter in terms of the expected return versus risk (variance) for a given asset data. The Optimize API has a method for obtaining efficient frontiers by successively calculting optimal portfolios for varying levels of return. By default, the smallest return is the one attained by the minimum variance portfolio, and the largest return is that of the largest individual asset.
# 
# Calling the frontier api method is very simlar to performing a single portfolio optimization. The following call:

# In[45]:


# vector of expected returns
r = assets['Return (%)'].values

# covariance matrix
sigmas = np.diag(assets['STD (%)'].values)
Q = sigmas @ cov.values @ sigmas

# options
options = {
    'number_of_points': 20
}

# prepare data for submitting request to api
data = {
    'Q': Q.tolist(),
    'r': r.tolist(),
    'options': options
}
sol = client.call('frontier', data)


# will produce 20 points on the efficient frontier of the classic Markowitz portfolio optimization propblem, that is one without any additional constraints on the portfolio or the positions:

# In[46]:


frontier = sol['frontier']
m = len(frontier)
mu0 = np.zeros((m,))
var0 = np.zeros((m,))
pos0 = np.zeros((m, n))
for i, e in enumerate(frontier):
    point = e['sol']
    mu0[i] = e['mu']
    var0[i] = point['obj']
    pos0[i, :] = point['x']
    if point['obj'] is not None:
        print('  {}. mu = {:.3f}%, std = {:.3f}% ({})'.format(i+1, mu0[i], math.sqrt(var0[i]), point['status']))
    else:
        print('  {}. mu = {:.3f}%, std = --- ({})'.format(i+1, mu0[i], point['status']))


# These points can be combine to produce a plot of the frontier:

# In[47]:


rmax = np.max(r)
rmin = np.min((np.min(r), np.min(mu0)))
dr = rmax - rmin
colors = np.array(plt.rcParams['axes.prop_cycle'].by_key()['color'])


# In[48]:


# plot frontier
xlim = ((rmin-0.1*dr), (rmax+0.1*dr))

for i, label in enumerate(assets['Label'].values): 
    plt.plot(np.sqrt(Q[i,i]), r[i], 'x', color=colors[np.mod(i, len(colors))])
plt.plot(np.sqrt(var0), mu0, 'k')
plt.xlabel('std (%)')
plt.ylabel('return (%)')
plt.title('Efficient frontier')
plt.ylim(xlim)
plt.legend(assets['Label'].values)
plt.grid()
plt.show()


# ## The frontier without short sales
# 
# The Optimize API let's one add any constraint to the calculation of the points in the frontier.
# 
# For example:

# In[49]:


# vector of expected returns
r = assets['Return (%)'].values

# covariance matrix
sigmas = np.diag(assets['STD (%)'].values)
Q = sigmas @ cov.values @ sigmas

# options
options = {
    'number_of_points': 20,
    'short': False
}

# prepare data for submitting request to api
data = {
    'Q': Q.tolist(),
    'r': r.tolist(),
    'options': options
}
sol = client.call('frontier', data)


# will produce 20 points on the efficient frontier obtained for portfolios for which short sales are not permitted: 

# In[50]:


frontier = sol['frontier']
m = len(frontier)
mu1 = np.zeros((m,))
var1 = np.zeros((m,))
pos1 = np.zeros((m, n))
for i, e in enumerate(frontier):
    point = e['sol']
    mu1[i] = e['mu']
    var1[i] = point['obj']
    pos1[i, :] = point['x']
    if point['obj'] is not None:
        print('  {}. mu = {:.3f}%, std = {:.3f}% ({})'.format(i+1, mu1[i], math.sqrt(var1[i]), point['status']))
    else:
        print('  {}. mu = {:.3f}%, std = --- ({})'.format(i+1, mu1[i], point['status']))


# Again these points can be combined to produce a plot of the frontier:

# In[51]:


rmax = np.max(r)
rmin = np.min((np.min(r), np.min(mu0)))
dr = rmax - rmin
colors = np.array(plt.rcParams['axes.prop_cycle'].by_key()['color'])


# In[52]:


# plot frontier
ylim = ((rmin-0.1*dr), (rmax+0.1*dr))

for i, label in enumerate(assets['Label'].values): 
    plt.plot(np.sqrt(Q[i,i]), r[i], 'x', color=colors[np.mod(i, len(colors))])
plt.plot(np.sqrt(var0), mu0, 'k--')
plt.plot(np.sqrt(var1), mu1, 'k')
plt.xlabel('std (%)')
plt.ylabel('return (%)')
plt.title('Efficient frontier')
plt.ylim(ylim)
plt.legend(assets['Label'].values)
plt.grid()
plt.show()


# We have combined in this plot also the standard frontier with short sales (dashed) for comparison.

# ## The frontier with general constraints

# All options available for calculation of optimal portfolios also apply to the frontier.
# 
# For example, the requirement of no short sales plus at least 50% on S & P and Bonds can be applied to the frontier as follows:

# In[53]:


# options
options = {
    'number_of_points': 20,
    'short': False 
}

# prepare data for submitting request to api
data = {
    'Q': Q.tolist(),
    'r': r.tolist(),
    'mu': mu,
    'labels': assets['Label'].tolist(),
    'constraints': [{'label': 's&p_bonds', 'set': ['S & P', 'Bonds'], 'bounds': [.5, np.inf]}],
    'options': options
}
sol = client.call('frontier', data)


# In[54]:


frontier = sol['frontier']
m = len(frontier)
mu2 = np.zeros((m,))
var2 = np.zeros((m,))
pos2 = np.zeros((m, n))
for i, e in enumerate(frontier):
    point = e['sol']
    mu2[i] = e['mu']
    var2[i] = point['obj']
    pos2[i, :] = point['x']
    if point['obj'] is not None:
        print('  {}. mu = {:.3f}%, std = {:.3f}% ({})'.format(i+1, mu2[i], math.sqrt(var2[i]), point['status']))
    else:
        print('  {}. mu = {:.3f}%, std = --- ({})'.format(i+1, mu2[i], point['status']))


# Note how the additional constraints lead to infeasible portfolios for higher levels of return.
# 
# This frontier's graphical representation can be calculated as in:

# In[55]:


rmax = np.max(r)
rmin = np.min((np.min(r), np.min(mu0)))
dr = rmax - rmin
colors = np.array(plt.rcParams['axes.prop_cycle'].by_key()['color'])


# In[56]:


# plot frontier
ylim = ((rmin-0.1*dr), (rmax+0.1*dr))

for i, label in enumerate(assets['Label'].values): 
    plt.plot(np.sqrt(Q[i,i]), r[i], 'x', color=colors[np.mod(i, len(colors))])
plt.plot(np.sqrt(var0), mu0, 'k--')
plt.plot(np.sqrt(var1), mu1, 'k:')
plt.plot(np.sqrt(var2), mu2, 'k')
plt.xlabel('std (%)')
plt.ylabel('return (%)')
plt.title('Efficient frontier')
plt.ylim(ylim)
plt.legend(assets['Label'].values)
plt.grid()
plt.show()


# ## A strongly constrained frontier
# 
# Consider now the cashflow problem we discussed above in which one would like to invest additional funds to improve their returns without selling any assets. The corresponding frontier can be obtained as follows:

# In[57]:


# initial portfolio
x0 = assets['Initial (%)'].values/100

# cashflow
cashflow = 20/100

# options
options = {
    'number_of_points': 20,
    'sell': False
}

# prepare data for submitting request to api
data = {
    'Q': Q.tolist(),
    'r': r.tolist(),
    'x0': x0.tolist(),
    'cashflow': cashflow,
    'options': options
}
sol = client.call('frontier', data)


# Note that in this case many returns lead to infeasible portfolios, and the efficient frontier covers a much smaller range of returns:

# In[58]:


frontier = sol['frontier']
m = len(frontier)
mu3 = np.zeros((m,))
var3 = np.zeros((m,))
pos3 = np.zeros((m, n))
for i, e in enumerate(frontier):
    point = e['sol']
    mu3[i] = e['mu']
    var3[i] = point['obj']
    pos3[i, :] = point['x']
    if point['obj'] is not None:
        print('  {}. mu = {:.3f}%, std = {:.3f}% ({})'.format(i+1, mu3[i], math.sqrt(var3[i]), point['status']))
    else:
        print('  {}. mu = {:.3f}%, std = --- ({})'.format(i+1, mu3[i], point['status']))


# as visualized in the following plot:

# In[59]:


rmax = np.max(r)
rmin = np.min((np.min(r), np.min(mu0)))
dr = rmax - rmin
colors = np.array(plt.rcParams['axes.prop_cycle'].by_key()['color'])


# In[60]:


# plot frontier
ylim = ((rmin-0.1*dr), (rmax+0.1*dr))

for i, label in enumerate(assets['Label'].values): 
    plt.plot(np.sqrt(Q[i,i]), r[i], 'x', color=colors[np.mod(i, len(colors))])
plt.plot(np.sqrt(var0), mu0, 'k--')
plt.plot(np.sqrt(var1), mu1, 'k:')
plt.plot(np.sqrt(var2), mu2, 'k-.')
plt.plot(np.sqrt(var3), mu3, 'k')
plt.xlabel('std (%)')
plt.ylabel('return (%)')
plt.title('Efficient frontier')
plt.ylim(ylim)
plt.legend(assets['Label'].values)
plt.grid()
plt.show()


# Zooming in:

# In[61]:


# plot frontier
ylim = (11.5, 13.5)
xlim = (9, 13)

for i, label in enumerate(assets['Label'].values): 
    plt.plot(np.sqrt(Q[i,i]), r[i], 'x', color=colors[np.mod(i, len(colors))])
plt.plot(np.sqrt(var0), mu0, 'k--')
plt.plot(np.sqrt(var1), mu1, 'k:')
plt.plot(np.sqrt(var2), mu2, 'k-.')
plt.plot(np.sqrt(var3), mu3, 'k')
plt.xlabel('std (%)')
plt.ylabel('return (%)')
plt.title('Efficient frontier')
plt.ylim(ylim)
plt.xlim(xlim)
plt.grid()
plt.show()


# helps to visualize the three frontiers.
# 
# Another interesting visualizion is that of the resulting trades. The following plot shows the frontier, the corresponding positions as a percentage of the total holdings and the corresponding trades as a percentage of the cashflow.

# In[62]:


mu = mu3
var = var3
pos = pos3
trades = (pos - x0) / cashflow

# plot frontier
plt.figure(figsize=(10, 10))

ax = plt.subplot(3, 1, 1)
for i, label in enumerate(assets['Label'].values): 
    plt.plot(r[i], np.sqrt(Q[i,i]), 'x', color=colors[np.mod(i, len(colors))])
ax.plot(mu0, np.sqrt(var0), 'k--')
ax.plot(mu1, np.sqrt(var1), 'k:')
ax.plot(mu2, np.sqrt(var2), 'k-.')
ax.plot(mu, np.sqrt(var), 'k')
ax.set_ylabel('std (%)')
ax.set_title('Efficient frontier')
ax.set_xlim(ylim)
ax.set_ylim(xlim)
ax.grid()

ax = plt.subplot(3, 1, 2)
ax.plot(mu, 100*pos)
ax.set_ylabel('positions')
ax.set_title('Positions')
ax.set_xlim(ylim)
ax.legend(assets['Label'].values)
ax.grid()

ax = plt.subplot(3, 1, 3)
ax.plot(mu, 100*trades)
ax.set_xlabel('return (%)')
ax.set_ylabel('trades')
ax.set_title('Buy/Sell')
ax.set_xlim(ylim)
ax.grid()

plt.subplots_adjust(hspace=.4)
plt.show()


# Note that the frontier in this case is represented with the coordinates reversed, that is risk versus return, as opposed to the standard return versus risk.