import numpy
from cvxopt.base import spmatrix
from cvxopt.base import matrix as _matrix
from cvxopt.blas import dot
from cvxopt.solvers import qp, options
RealNumber = float
Integer=int
|
|
n=4
n_S=numpy.matrix([[4e-2, 6e-3, -4e-3, 0.0],[6e-3, 1e-2, 0.0, 0.0],[-4e-3, 0.0, 2.5e-3, 0.0],[0.0, 0.0, 0.0, 0.0]])
|
|
S=_matrix(n_S)
|
|
n_pbar=numpy.matrix([0.12,0.10,0.07,0.03])
n_pbar.shape=(4,1)
pbar = _matrix(n_pbar)
|
|
G=_matrix(0.0,(n,n))
G[::n+1]=-1.0
#print(G)
|
|
h = _matrix(0.0, (n,1))
A = _matrix(1.0, (1,n))
b = _matrix(1.0)
|
|
N = 100
mus = [ 10.0^(5.0*t/N-1.0) for t in xrange(N) ]
options['show_progress'] = False
xs = [qp(float(mu)*S, -pbar, G, h, A, b)['x'] for mu in mus ]
|
|
returns = [ dot(pbar,x) for x in xs ]
risks = [ sqrt(dot(x, S*x)) for x in xs ]
|
|
P1=list_plot(zip(risks,returns), plotjoined=True)
P1.axes_labels(['standard deviation','expected return'])
P1.show()
|
|
c1 = [ x[0] for x in xs ]
c2 = [ x[0] + x[1] for x in xs ]
c3 = [ x[0] + x[1] + x[2] for x in xs ]
c4 = [ x[0] + x[1] + x[2] + x[3] for x in xs ]
|
|
import pylab
pylab.figure(1)
pylab.fill(risks + [.20], c1 + [0.0], facecolor = '#F0F0F0')
pylab.fill(risks[-1::-1] + risks, c2[-1::-1] + c1, facecolor = '#D0D0D0')
pylab.fill(risks[-1::-1] + risks, c3[-1::-1] + c2, facecolor = '#F0F0F0')
pylab.fill(risks[-1::-1] + risks, c4[-1::-1] + c3, facecolor = '#D0D0D0')
pylab.axis([0.0, 0.2, 0.0, 1.0])
pylab.savefig('temp.png')
|
from cvxopt import blas, lapack
from cvxopt.solvers import cp
def spdiag(x):
return spmatrix(x,range(len(x)),range(len(x)))
my_list=[1.0,1.0,1.0,1.0]
print(len(my_list))
spdiag(my_list)
4 <4x4 spmatrix, tc='d', nnz=4> 4 <4x4 spmatrix, tc='d', nnz=4> |
V = _matrix([-2.1213, 2.1213,
-2.2981, 1.9284,
-2.4575, 1.7207,
-2.5981, 1.5000,
-2.7189, 1.2679,
-2.8191, 1.0261,
-2.8978, 0.7765,
-2.9544, 0.5209,
-2.9886, 0.2615,
-3.0000, 0.0000,
1.5000, 0.0000,
1.4772, -0.2605,
1.4095, -0.5130,
1.2990, -0.7500,
1.1491, -0.9642,
0.9642, -1.1491,
0.7500, -1.2990,
0.5130, -1.4095,
0.2605, -1.4772,
0.0000, -1.5000 ], (2,20))
|
|
n = V.size[1]
G = spmatrix(-1.0, range(n), range(n))
h = _matrix(0.0, (n,1))
A = _matrix(1.0, (1,n))
b = _matrix(1.0)
|
|
# D-design
#
# minimize f(x) = -log det V*diag(x)*V'
# subject to x >= 0
# sum(x) = 1
#
# The gradient and Hessian of f are
#
# gradf = -diag(V' * X^-1 * V)
# H = (V' * X^-1 * V)**2.
#
# where X = V * diag(x) * V'.
def F(x=None, z=None):
if x is None:
return 0, _matrix(float(1.0), (n,1))
#print(x)
#print(spdiag(x))
X = V * spdiag(x) * V.T
L = +X
try:
lapack.potrf(L)
except ArithmeticError:
return None
f = -2.0 * (log(L[0,0]) + log(L[1,1]))
W = +V
blas.trsm(L, W)
gradf = _matrix(-1.0, (1,2)) * W**2
if z is None:
return f, gradf
H = _matrix(0.0,(n,n))
blas.syrk(W, H, trans='T')
return f, gradf, z[0] * H**2
|
|
xd = cp(F, G, h, A = A, b = b)['x']
|
|
pylab.figure(2, facecolor='w', figsize=(6,6))
pylab.plot(V[0,:], V[1,:],'ow', [0], [0], 'k+')
I = [ k for k in xrange(n) if xd[k] > 1e-5 ]
pylab.plot(V[0,I], V[1,I],'or')
[<matplotlib.lines.Line2D object at 0x4c55450>, <matplotlib.lines.Line2D object at 0x4c55490>] [<matplotlib.lines.Line2D object at 0x4c55450>, <matplotlib.lines.Line2D object at 0x4c55490>] |
from numpy import array as _array
# Enclosing ellipse is {x | x' * (V*diag(xe)*V')^-1 * x = sqrt(2)}
nopts = 1000
angles = _array(_matrix([ float(c) for c in [ a*2.0*pi/nopts for a in xrange(nopts)]] , (1,nopts) ))
circle = _matrix(0.0, (2,nopts))
circle[0,:], circle[1,:] = cos(angles), sin(angles)
|
|
W = V * spdiag(xd) * V.T
lapack.potrf(W)
ellipse = sqrt(2.0) * circle
blas.trmm(W, ellipse)
pylab.plot(ellipse[0,:].T, ellipse[1,:].T, 'k--')
pylab.axis([-5, 5, -5, 5])
pylab.title('D-optimal design (fig. 7.9)')
pylab.axis('off')
pylab.savefig('temp2.png')
|
|
|
|
|
|
|
|
|