Clase 11: Algunas mejoras a los códigos para simular y optimizar portafolios

Juan Diego Sánchez Torres,

Profesor, MAF ITESO

  • Departamento de Matemáticas y Física
  • dsanchez@iteso.mx
  • Tel. 3669-34-34 Ext. 3069
  • Oficina: Cubículo 4, Edificio J, 2do piso

1. Motivación

En primer lugar, para poder bajar precios y información sobre opciones de Yahoo, es necesario cargar algunos paquetes de Python. En este caso, el paquete principal será Pandas. También, se usarán el Scipy y el Numpy para las matemáticas necesarias y, el Matplotlib y el Seaborn para hacer gráficos de las series de datos. Finalmente, se usará el paquete cvxopt para optimización convexa, para instalar ingrese en terminal la instrucción: conda install -c anaconda cvxopt


In [1]:
#importar los paquetes que se van a usar
import pandas as pd
import numpy as np
import datetime
from datetime import datetime
import scipy.stats as stats
import scipy as sp
import matplotlib.pyplot as plt
import seaborn as sns
import sklearn.covariance as skcov
import cvxopt as opt
from cvxopt import blas, solvers
solvers.options['show_progress'] = False
%matplotlib inline
pd.set_option('display.notebook_repr_html', True)
pd.set_option('display.max_columns', 6)
pd.set_option('display.max_rows', 10)
pd.set_option('display.width', 78)
pd.set_option('precision', 3)
#Funciones para portafolios
import portfolio_func
from pyomo.environ import *
infinity = float('inf')
import statsmodels.api as sm


C:\Users\DSANCHEZ\AppData\Local\Continuum\Anaconda3\lib\site-packages\statsmodels\compat\pandas.py:56: FutureWarning: The pandas.core.datetools module is deprecated and will be removed in a future version. Please use the pandas.tseries module instead.
  from pandas.core import datetools

In [2]:
assets = ['AAPL','MSFT','AA','AMZN','KO','QAI']
closes = portfolio_func.get_historical_closes(assets, '2016-01-01', '2017-09-22')

In [3]:
daily_returns=portfolio_func.calc_daily_returns(closes)
huber = sm.robust.scale.Huber()
#Mean and standar deviation returns
returns_av, scale = huber(daily_returns)

In [21]:
model.assets = Set()
model.T = Set(initialize = range(1994, 2014))
model.max_risk = Param(mutable = True, initialize = .00305)
model.R = Param(model.T, model.assets)

def mean_init(model, j):
    return sum(model.R[i, j] for i in model.T)/len(model.T)
model.mean = Param(model.assets, initialize = mean_init)

def Q_init(model, i, j):
    return sum((model.R[k, i] - model.mean[i])*(model.R[k, j] - model.mean[j]) for k in model.T)
model.Q = Param(model.assets, model.assets, initialize = Q_init)
model.alloc = Var(model.assets, within=NonNegativeReals)

def risk_bound_rule(model):
    return (sum(sum(model.Q[i, j] * model.alloc[i] * model.alloc[j] for i in model.assets)for j in model.assets) <= model.max_risk)
model.risk_bound = Constraint(rule=risk_bound_rule)

def tot_mass_rule(model):
    return (sum(model.alloc[j] for j in model.assets) == 1)
model.tot_mass = Constraint(rule=tot_mass_rule)

def objective_rule(model):
    return summation(model.mean, model.alloc)
model.objective = Objective(sense=maximize, rule=objective_rule)


WARNING: Implicitly replacing the Component attribute assets (type=<class 'pyomo.core.base.param.IndexedParam'>) on block unknown with a new Component (type=<class 'pyomo.core.base.sets.SimpleSet'>).
	This is usually indicative of a modelling error.
	To avoid this warning, use block.del_component() and block.add_component().
WARNING: Implicitly replacing the Component attribute T (type=<class 'pyomo.core.base.var.IndexedVar'>) on block unknown with a new Component (type=<class 'pyomo.core.base.sets.SimpleSet'>).
	This is usually indicative of a modelling error.
	To avoid this warning, use block.del_component() and block.add_component().
WARNING: Implicitly replacing the Component attribute mean (type=<class 'pyomo.core.base.param.IndexedParam'>) on block unknown with a new Component (type=<class 'pyomo.core.base.param.IndexedParam'>).
	This is usually indicative of a modelling error.
	To avoid this warning, use block.del_component() and block.add_component().
ERROR: Constructing component 'tot_mass' from data=None failed:
	ValueError: Invalid constraint expression. The constraint expression resolved to a trivial Boolean (False) instead of a Pyomo object. Please modify your rule to return Constraint.Infeasible instead of False.
	
	Error thrown for Constraint 'tot_mass'
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-21-0ba2a64614e3> in <module>()
     19 def tot_mass_rule(model):
     20     return (sum(model.alloc[j] for j in model.assets) == 1)
---> 21 model.tot_mass = Constraint(rule=tot_mass_rule)
     22 
     23 def objective_rule(model):

~\AppData\Local\Continuum\Anaconda3\lib\site-packages\pyomo\core\base\block.py in __setattr__(self, name, val)
    540                 # Pyomo components are added with the add_component method.
    541                 #
--> 542                 self.add_component(name, val)
    543             else:
    544                 #

~\AppData\Local\Continuum\Anaconda3\lib\site-packages\pyomo\core\base\block.py in add_component(self, name, val)
    967                               _blockName, str(data) )
    968             try:
--> 969                 val.construct(data)
    970             except:
    971                 err = sys.exc_info()[1]

~\AppData\Local\Continuum\Anaconda3\lib\site-packages\pyomo\core\base\constraint.py in construct(self, data)
    736 
    737             assert None not in self._data
--> 738             cdata = self._check_skip_add(None, tmp, condata=self)
    739             if cdata is not None:
    740                 # this happens as a side-effect of set_value on

~\AppData\Local\Continuum\Anaconda3\lib\site-packages\pyomo\core\base\constraint.py in _check_skip_add(self, index, expr, condata)
    902                        expr and "Feasible" or "Infeasible",
    903                        expr,
--> 904                        self._data[index].name))
    905 
    906         #

ValueError: Invalid constraint expression. The constraint expression resolved to a trivial Boolean (False) instead of a Pyomo object. Please modify your rule to return Constraint.Infeasible instead of False.

Error thrown for Constraint 'tot_mass'

In [5]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [6]:
!type dietdata.dat


param:  F:                          c     V  :=
  "Cheeseburger"                 1.84   4.0  
  "Ham Sandwich"                 2.19   7.5  
  "Hamburger"                    1.84   3.5  
  "Fish Sandwich"                1.44   5.0  
  "Chicken Sandwich"             2.29   7.3  
  "Fries"                         .77   2.6  
  "Sausage Biscuit"              1.29   4.1  
  "Lowfat Milk"                   .60   8.0 
  "Orange Juice"                  .72  12.0 ;

param Vmax := 75.0;

param:  N:       Nmin   Nmax :=
        Cal      2000      .
        Carbo     350    375
        Protein    55      .
        VitA      100      .
        VitC      100      .
        Calc      100      .
        Iron      100      . ;

param a:
                               Cal  Carbo Protein   VitA   VitC  Calc  Iron :=
  "Cheeseburger"               510     34     28     15      6    30    20
  "Ham Sandwich"               370     35     24     15     10    20    20
  "Hamburger"                  500     42     25      6      2    25    20
  "Fish Sandwich"              370     38     14      2      0    15    10
  "Chicken Sandwich"           400     42     31      8     15    15     8
  "Fries"                      220     26      3      0     15     0     2
  "Sausage Biscuit"            345     27     15      4      0    20    15
  "Lowfat Milk"                110     12      9     10      4    30     0
  "Orange Juice"                80     20      1      2    120     2     2 ;

In [7]:
!pyomo solve --solver=glpk diet.py dietdata.dat


[    0.00] Setting up Pyomo environment
[    0.00] Applying Pyomo preprocessing actions
[    0.01] Creating model
[    0.05] Applying solver
[    0.14] Processing results
    Number of solutions: 1
    Solution Information
      Gap: 0.0
      Status: optimal
      Function Value: 15.05
    Solver results file: results.yml
[    0.14] Applying Pyomo postprocessing actions
[    0.14] Pyomo Finished

In [8]:
!type results.yml


# ==========================================================
# = Solver Results                                         =
# ==========================================================
# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Name: unknown
  Lower bound: 15.05
  Upper bound: 15.05
  Number of objectives: 1
  Number of constraints: 10
  Number of variables: 10
  Number of nonzeros: 77
  Sense: minimize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  Termination condition: optimal
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: 23
      Number of created subproblems: 23
  Error rc: 0
  Time: 0.06200003623962402
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 1
  number of solutions displayed: 1
- Gap: 0.0
  Status: optimal
  Message: None
  Objective:
    cost:
      Value: 15.05
  Variable:
    x[Cheeseburger]:
      Value: 4
    x[Fish Sandwich]:
      Value: 1
    x[Fries]:
      Value: 5
    x[Lowfat Milk]:
      Value: 4
  Constraint: No values

In [ ]:

2. Uso de Pandas para descargar datos de precios de cierre

Una vez cargados los paquetes, es necesario definir los tickers de las acciones que se usarán, la fuente de descarga (Yahoo en este caso, pero también se puede desde Google) y las fechas de interés. Con esto, la función DataReader del paquete pandas_datareader bajará los precios solicitados.

Nota: Usualmente, las distribuciones de Python no cuentan, por defecto, con el paquete pandas_datareader. Por lo que será necesario instalarlo aparte. El siguiente comando instala el paquete en Anaconda: conda install -c conda-forge pandas-datareader


In [ ]:

3. Formulación del riesgo de un portafolio y simulación Montecarlo


In [ ]:


In [9]:
r=0.0001
results_frame = portfolio_func.sim_mont_portfolio(daily_returns,100000,r)

In [10]:
#Sharpe Ratio
max_sharpe_port = results_frame.iloc[results_frame['Sharpe'].idxmax()]
#Menor SD
min_vol_port = results_frame.iloc[results_frame['SD'].idxmin()]

In [11]:
plt.scatter(results_frame.SD,results_frame.Returns,c=results_frame.Sharpe,cmap='RdYlBu')
plt.xlabel('Volatility')
plt.ylabel('Returns')
plt.colorbar()
#Sharpe Ratio
plt.scatter(max_sharpe_port[1],max_sharpe_port[0],marker=(5,1,0),color='r',s=1000);
#Menor SD
plt.scatter(min_vol_port[1],min_vol_port[0],marker=(5,1,0),color='g',s=1000);



In [12]:
pd.DataFrame(max_sharpe_port)


Out[12]:
10881
Returns 0.210
SD 0.113
Sharpe 1.861
AA 0.048
AAPL 0.215
AMZN 0.150
KO 0.381
MSFT 0.077
QAI 0.130

In [13]:
pd.DataFrame(min_vol_port)


Out[13]:
89245
Returns 0.115
SD 0.076
Sharpe 1.501
AA 0.012
AAPL 0.089
AMZN 0.043
KO 0.243
MSFT 0.009
QAI 0.604

4. Optimización de portafolios


In [14]:
N=5000
results_frame_optim = portfolio_func.optimal_portfolio(daily_returns,N,r)

In [15]:
#Montecarlo
plt.scatter(results_frame.SD,results_frame.Returns,c=results_frame.Sharpe,cmap='RdYlBu')
plt.xlabel('Volatility')
plt.ylabel('Returns')
plt.colorbar()
#Markowitz
plt.plot(results_frame_optim.SD, results_frame_optim.Returns, 'b-o');



In [16]:
#Sharpe Ratio
max_sharpe_port_optim = results_frame_optim.iloc[results_frame_optim['Sharpe'].idxmax()]
#Menor SD
min_vol_port_optim = results_frame_optim.iloc[results_frame_optim['SD'].idxmin()]

In [17]:
#Markowitz
plt.scatter(results_frame_optim.SD,results_frame_optim.Returns,c=results_frame_optim.Sharpe,cmap='RdYlBu');
plt.xlabel('Volatility')
plt.ylabel('Returns')
plt.colorbar()
#Sharpe Ratio
plt.scatter(max_sharpe_port_optim[1],max_sharpe_port_optim[0],marker=(5,1,0),color='r',s=1000);
#SD
plt.scatter(min_vol_port_optim[1],min_vol_port_optim[0],marker=(5,1,0),color='g',s=1000);



In [18]:
pd.DataFrame(max_sharpe_port_optim)


Out[18]:
2216
Returns 0.211
SD 0.113
Sharpe 1.862
AA 0.049
AAPL 0.234
AMZN 0.149
KO 0.354
MSFT 0.078
QAI 0.138

In [19]:
pd.DataFrame(min_vol_port_optim)


Out[19]:
4999
Returns 0.091
SD 0.074
Sharpe 1.222
AA 0.002
AAPL 0.043
AMZN 0.007
KO 0.217
MSFT 0.021
QAI 0.711