%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)
## NUMBER OF ASSETS
n_assets = 4
## NUMBER OF OBSERVATIONS
n_obs = 1000
return_vec = np.random.randn(n_assets, n_obs) + 0.000001
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
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)
0.0555080150868
<matplotlib.legend.Legend at 0x1102b0450>