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 [ ]: