Notebook
In [1]:
import numpy as np
import matplotlib.pyplot as plt
import cvxopt as opt
from cvxopt import blas, solvers
import pandas as pd

np.random.seed(123)

# Turn off progress printing 
solvers.options['show_progress'] = False
In [2]:
# get pricing data where symbols = ['xxx','xxx'],dates = 'yyyy-mm-dd'
def get_stock_pricing (symbols,startdate, enddate):
    return get_pricing(symbols,
                     start_date=startdate, end_date=enddate)
In [3]:
# get returns
def get_returns(symbols,startdate,enddate):
    prices = get_pricing(symbols,fields='close_price',start_date=startdate, end_date=enddate,frequency='daily')
    returns = prices.pct_change().dropna()*100
    return_vec=np.array(returns)
    return return_vec.T
In [ ]:
 
In [101]:
#########################################################################
#                     HISTORICAL SERIES INPUT
#########################################################################
symbols = ['SPY','IWD','VNQ','IWM','HEDJ']

#['AAPL','TSLA','JPM','MSFT']

startdate = '2008-6-30'
enddate = '2016-9-30'
#get_returns(symbols,startdate,enddate)    
In [102]:
# define random weights where n here is number of assets
def rand_weights(n):
    ''' Produces n random weights that sum to 1 '''
    k = np.random.rand(n)
    return k / sum(k)
In [103]:
# define random portfolios of the randomly weighted given assets
def random_portfolio(returns):
    ''' 
    Returns the mean and standard deviation of returns for a random portfolio
    '''

    p = np.asmatrix(np.mean(returns, axis=1))
    w = np.asmatrix(rand_weights(returns.shape[0]))
    C = np.asmatrix(np.cov(returns))
    
    mu = w * p.T
    sigma = np.sqrt(w * C * w.T)
    
    # This recursion reduces outliers to keep plots pretty
    if sigma > 3:
        return random_portfolio(returns)
    return mu, sigma
In [104]:
return_vec = get_returns(symbols,startdate,enddate)

n_portfolios = 500

means, stds = np.column_stack([random_portfolio(return_vec) for _ in xrange(n_portfolios)])
In [105]:
plt.plot(stds, means, 'o', markersize=5)
plt.xlabel('std')
plt.ylabel('mean')
plt.title('Mean and standard deviation of returns of randomly generated portfolios');
In [106]:
def optimal_portfolio(returns, mus):
    n = len(returns)
    returns = np.asmatrix(returns)
    
    #mus = [0.00, 0.005, 0.010, 0.015, 0.020, 0.025, 0.030, 0.040, 0.050, 0.060, 0.070, 0.1]
    
    # Convert to cvxopt matrices
    S = opt.matrix(np.cov(returns))
    pbar = opt.matrix(np.mean(returns, axis=1))
    
    # Create constraint matrices
    G = -opt.matrix(np.eye(n))   # negative n x n identity matrix
    h = opt.matrix(0.0, (n ,1))
    A = opt.matrix(1.0, (1, n))
    b = opt.matrix(1.0)
    
    # Calculate efficient frontier weights using quadratic programming
    portfolios = [solvers.qp(mu*S, -pbar, G, h, A, b)['x'] 
                  for mu in mus]
        
    ## CALCULATE RISKS AND RETURNS FOR FRONTIER
    returns = [blas.dot(pbar, x) for x in portfolios]
    risks = [np.sqrt(blas.dot(x, S*x)) for x in portfolios]
    ## CALCULATE THE 2ND DEGREE POLYNOMIAL OF THE FRONTIER CURVE
    m1 = np.polyfit(returns, risks, 2)
    x1 = np.sqrt(m1[2] / m1[0])
    # CALCULATE THE OPTIMAL PORTFOLIO
    wt = solvers.qp(opt.matrix(x1 * S), -pbar, G, h, A, b)['x'] 
      
    return np.asarray(wt), returns, risks
In [107]:
    
N = 100
mus = [10**(5.0 * t/N - 1.0) for t in range(N)]

weights, returns, risks = optimal_portfolio(return_vec, mus)

print weights

plt.plot(stds, means, 'o')
plt.ylabel('mean')
plt.xlabel('std')
plt.plot(risks, returns, 'y-o');
[[  8.00354698e-01]
 [  3.30577373e-08]
 [  1.99644641e-01]
 [  6.23524810e-07]
 [  4.72954156e-09]]
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]: