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

# In[5]:


get_ipython().run_line_magic('matplotlib', 'inline')
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = (10.0, 7.0)
from scipy.optimize import curve_fit
import cvxopt as opt
from cvxopt import blas, solvers
solvers.options['show_progress'] = False
np.random.seed(123)


# In[6]:


## NUMBER OF ASSETS
n_assets = 4

## NUMBER OF OBSERVATIONS
n_obs = 1000

return_vec = np.random.randn(n_assets, n_obs) + 0.000001


# In[7]:


def rand_weights(n):
    ''' Produces n random weights that sum to 1 '''
    k = np.random.rand(n)
    return k / sum(k)

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 > 2:
        return random_portfolio(returns)
    return mu, sigma

def optimal_portfolio(returns):
    n = len(returns)
    returns = np.asmatrix(returns)
    
    N = 100
    mus = [10**(5.0 * t/N - 1.0) for t in range(N)]
    
    # 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])
    print x1
    # CALCULATE THE OPTIMAL PORTFOLIO
    wt = solvers.qp(opt.matrix(x1 * S), -pbar, G, h, A, b)['x']
    return np.asarray(wt), returns, risks



# In[8]:


n_portfolios = 500
means, stds = np.column_stack([
    random_portfolio(return_vec) 
    for _ in xrange(n_portfolios)
])

weights, returns, risks = optimal_portfolio(return_vec)

# Quadratic Fit
quadratic = np.poly1d(np.polyfit(returns, risks, 2))

# Try fitting a Hyperbola
def hyperbola(x, x0, y0, a, b):
    """
    formula used: 
        (y - y0)^2 / a^2 - (x - x0)^2 / b^2 = 1
    """
    return y0**2 + np.sqrt(a**2 * ((x - x0)**2 / b**2 + 1))


# Initial params guess, these seem pretty sensitive
guess = 0.05, 0.05, np.mean(risks) / 2, np.mean(returns) / 2

params, _ = curve_fit(hyperbola, returns, risks, guess)

plt.plot(means, stds, 'r+')
plt.xlabel('mean')
plt.ylabel('std')
plt.plot(returns, risks,  'y-o', label='Efficient Frontier')

plt.plot(returns, quadratic(returns), 'b-+', label='quadratic')
plt.plot(returns, hyperbola(returns, *params), 'g-x', label='hyperbolic')
plt.legend(loc=0)


# In[ ]: