Notebook
In [24]:
import math
import numpy as np
from cvxopt import matrix, spmatrix
from cvxopt.blas import dot
from cvxopt.solvers import sdp, conelp
from cvxopt import solvers
import scipy.stats as st
solvers.options['maxiters'] = 1000
solvers.options['show_progress'] = False

stocks = symbols(['SPY','VXX','EEM','GDX','XLF','UWTI','EWJ','IWM','QQQ','USO','XIV','TVIX','UVXY','UGAZ','FXI',
'EFA','XLE','VWO','EWZ','AMLP','XOP','XLV','XLI','JNK'])#,'HYG','SDS','XLU','XLK','RSX','DUST',
#'XLP','IYR','VEA','EWT','TLT','GDXJ','NUGT','UCO','XLY','EZU'])

prices = get_pricing(stocks,start_date='2013-01-01',end_date = '2013-06-01',fields='price',frequency='daily')
prices = prices.dropna(axis=1)

In [25]:
def getH(n, k):
H = np.ones((3*n*n+1, 1))
return H

def getG(n, k):
# we have n x n variables in Y
# we introduct n x n new variables for L1 norm constraint
# we have one more variable for z >= 0

G = np.zeros((3*n*n+1, 2*n*n+1))

# contraints for > 0 for auxillary variables
for i in range(0, n*n):
G[i, n*n+i] = -1.0

# constraints for inequalities
for i in range(0,n*n):
G[n*n+i,i] = 1.0
G[n*n+i,n*n+i] = -1.0

# constraints for inequalities
for i in range(0, n*n):
G[2*n*n+i,i] = -1.0
G[2*n*n+i,n*n+i] = -1.0

# constraint for z > 0
G[-1, -1] = -1.0

return G

def getGs(n):
Gs = np.zeros((n*n, 2*n*n+1))
for i in range(0, n*n):
Gs[i,i] = -1.0
return Gs

def getHs(n):
Hs = np.ones((n,n)) * -.000000000001
return Hs

In [26]:
K = 15
pval = prices.values
(row, col) = np.shape(pval)
data = pval - np.mean(pval, axis=0)
Y = data[1:, :]
Ym = data[:-1, :]

Als = np.linalg.solve(np.dot(Ym.T, Ym), np.dot(Ym.T, Y))
B = np.dot(Y.T, Y)
M = np.dot(Als.T, Als)
(n,n) = np.shape(M)
c = np.ravel(M.T).tolist() # constants for Y i.e. n x n
c = c + ([0] * (n*n))      # constants for new variables W
c.append(0)                # constant for z
c = np.asarray(c)

a1 = np.ravel(np.eye(n)).tolist() # for Y
a1 = a1 + ([0] * (n*n))           # for W
a1.append(-1.0)                   # for z

a2 = np.ravel(B.T).tolist()       # for Y
a2 = a2 + ([0] * (n*n+1))         # for W and z

a3 = ([0] * (n*n))                # for Y
a3 = a3 + ([1.0] * (n*n))         # for W
a3.append(-K)                     # for z

a4 = ([0] * (n*n))                # for Y
a4 = a4 + ([1.0] * (n*n))         # for W
a4.append(.0)                     # for z

Alist = [np.asarray(a1),np.asarray(a2), np.asarray(a3)]#, np.asarray(a4)]

Buff = ([0] * (n*n))
Buff = matrix(Buff, (n,n))

# condition for symmetry
for i in range(0, n):
for j in range(0, n):
if i <> j and Buff[i,j] == 0:
a = ([0] * (n*n))
m = matrix(a, (n,n))
Buff[i,j] = 1
Buff[j,i] = 1
m[i,j] = 1
m[j,i] = -1
a = np.ravel(m).tolist()
a = a + ([0] * (n*n))
a.append(0)
Alist.append(np.asarray(a))

A = np.matrix(Alist)

H = getH(n, K)
G = getG(n, K)
H1 = getHs(n)
G1 = getGs(n)

In [27]:
#CVXOPT begins
cm = matrix(c)
hm = matrix(H)
Gm = matrix(G)
G1 = [matrix(G1)]
H1 = [matrix(H1)]
Am = matrix(A)
bm = [0.0,1.0, 0]

Buff = ([0] * (n*n))
Buff = matrix(Buff, (n,n))

for i in range(0, n):
for j in range(0, n):
if i <> j and Buff[i,j] == 0:
bm.append(0)
Buff[i,j] = 1
Buff[j,i] = 1
bm = matrix(bm)
sol = sdp(c=cm,Gl=Gm,hl=hm,Gs=G1,hs=H1,A=Am,b=bm)
#sol = conelp(c=cm,G=Gm,h=hm,A=Am,b=bm)
print sol['status']

optimal

In [28]:
import scipy as sp
Y = matrix(sol['x'][0:(n*n)],(n,n))
z = sol['x'][-1]
evalue, evec =  sp.linalg.eig(Y)
print evalue

[  9.92080280e-01+0.j   9.80403479e-07+0.j   1.35541551e-08+0.j
8.63112281e-09+0.j   2.86141840e-09+0.j   2.26196256e-09+0.j
8.02364617e-10+0.j   2.17910461e-10+0.j   1.69918692e-10+0.j
8.46124257e-11+0.j   4.18403021e-11+0.j   3.95142285e-11+0.j
2.34291251e-11+0.j   1.77836397e-11+0.j   1.65797076e-11+0.j
1.38516740e-11+0.j   9.52635825e-12+0.j   7.08498646e-12+0.j
4.36129323e-12+0.j   2.45436254e-12+0.j   5.34756851e-13+0.j
-1.83169240e-14+0.j   2.52602039e-14+0.j]

In [29]:
idx = np.abs(evalue).argsort()[::-1]
print evalue[idx[0]]

(0.992080280103+0j)

In [30]:
W = np.ravel(evec[:, idx[0]])

In [31]:
print np.sort(np.abs(W))

[ 0.0013131   0.00952048  0.01375215  0.03241941  0.03735581  0.03892889
0.04480802  0.04964485  0.05803765  0.0837751   0.0971082   0.11789636
0.13245291  0.15778581  0.16155223  0.18027711  0.18494789  0.20743798
0.21357839  0.22876663  0.32325139  0.40801937  0.64060248]

In [32]:
import matplotlib.pyplot as plt
plt.plot(np.dot(pval, W * 1000))

Out[32]:
[<matplotlib.lines.Line2D at 0x7f261ae538d0>]
In [ ]: