In [25]:

    

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline


    

In [289]:

    

W = [1/3, 1/3, 1/3]
R = [.06, .02, .04]
C = np.matrix([[8, -2, 4], [-2, 2, -2], [4, -2, 8]])/1000
SD = [np.sqrt(C[i,j]) for i in range(len(R)) for j in range(len(R)) if i == j]
r = .05    # an initial target return
rf = .01   # risk free rate


    

In [290]:

    

def make_Muvec(R):        # R is a list [] of returns as real numbers
    Muvec = np.array(R)
    return Muvec
    
def make_Onevec(R):
    Onevec = np.ones(len(R))
    return Onevec

def make_Xvec(W):         # W is a list [] of weights as real numbers 
    Xvec = np.array(W)
    return Xvec


    

In [293]:

    

Equal_return = np.sum(Muvec * Xvec)
Equal_var = (np.matrix(Xvec) * C) * np.matrix(Xvec).T
Equal_vol = np.sqrt(Equal_var)
Equal_return, Equal_vol


    

    
Out[293]:





(0.039999999999999994, matrix([[ 0.04472136]]))

In [294]:

    

def make_LMtarget(r, R):      # r is a target return as a real number
    LMtarget = np.append(np.zeros(len(R)), [r,1])
    return np.matrix(LMtarget)


    

In [295]:

    

def Initialize(W, R, C, r, rf):
    M = make_Muvec(R)
    O = make_Onevec(R)
    X = make_Xvec(W)
    LMt = make_LMtarget(r,R)
    return M, O, X, LMt


    

In [296]:

    

Muvec, Onevec, Xvec, LMtarget = Initialize(W,R,C,r,rf)


    

In [297]:

    

def make_LMmat(C, R, Muvec, Onevec): # C is the covariance matrix as an NxN matrix of real numbers
    D = np.append(2*C, [Muvec], axis=0)
    E = np.append(D, [Onevec], axis=0)
    F = np.insert(E, len(R), np.append(Muvec, [0,0])*-1, axis=1)
    G = np.insert(F, len(R)+1, np.append(Onevec, [0,0])*-1, axis=1)
    return G


    

In [298]:

    

LMmat = make_LMmat(C, R, Muvec, Onevec)


    

In [299]:

    

def MinVar_target_return(C, LMmat, LMtarget):
    Xvec_target = LMmat.I * LMtarget.T
    Xvec_target = Xvec_target[:-2]
    target_var = Xvec_target.T * C * Xvec_target
    target_SD = np.sqrt(target_var)         # SD is standard deviation - measures volatility
    return Xvec_target, target_SD


    

In [300]:

    

Target_weights, target_SD = MinVar_target_return(C, LMmat, LMtarget)
Target_weights, target_SD


    

    
Out[300]:





(matrix([[ 0.76666667],
        [ 0.26666667],
        [-0.03333333]]),
 matrix([[ 0.06218253]]))

In [301]:

    

def MinVar_port(C, Muvec, Onevec):
    Muvec = np.matrix(Muvec).T
    Onevec = np.matrix(Onevec).T
    Xminvec = 1/np.sum(C.I * Onevec) * C.I * Onevec
    Xmin_return = 1/np.sum(C.I * Onevec) * Muvec.T * C.I * Onevec
    Xmin_SD = np.sqrt(Xminvec.T * C * Xminvec)
    return Xminvec, Xmin_return, Xmin_SD


    

In [302]:

    

Xminvec, Xmin_return, Xmin_SD = MinVar_port(C, Muvec, Onevec)
Xminvec, Xmin_return, Xmin_SD


    

    
Out[302]:





(matrix([[ 0.16666667],
        [ 0.66666667],
        [ 0.16666667]]),
 matrix([[ 0.03]]),
 matrix([[ 0.02581989]]))

In [303]:

    

def Sharpe_port(C, Muvec, rf):
    Muvec = np.matrix(Muvec).T
    Muhat = Muvec - rf
    XShvec = 1/np.sum(C.I * Muhat) * C.I * Muhat
    XSh_return = 1/np.sum(C.I * Muhat) * Muvec.T * C.I * Muhat
    XSh_SD = np.sqrt(1/np.sum(C.I * Muhat)**2 * Muhat.T * C.I * C * C.I * Muhat)
    Sharpe_ratio = (XSh_return - rf)/XSh_SD
    return XShvec, XSh_return, XSh_SD, Sharpe_ratio


    

In [304]:

    

XShvec, XSh_return, XSh_SD, Sharpe_ratio = Sharpe_port(C, Muvec, rf)
XShvec, XSh_return, XSh_SD, Sharpe_ratio


    

    
Out[304]:





(matrix([[ 0.29166667],
        [ 0.58333333],
        [ 0.125     ]]),
 matrix([[ 0.03416667]]),
 matrix([[ 0.02838231]]),
 matrix([[ 0.85146932]]))

In [305]:

    

def Find_Frontier(C, Muvec, rf):
    frontier = []
    frontier_weights = []
    r_points = np.arange(min(Muvec), max(Muvec)*1.1, .002)
    for r in r_points:
        LMtarget = make_LMtarget(r,R)
        Target_weights, target_SD = MinVar_target_return(C, LMmat, LMtarget)
        frontier.append([target_SD.item(), r])
        frontier_weights.append(Target_weights)
    return frontier, frontier_weights


    

In [306]:

    

frontier, frontier_weights = Find_Frontier(C,Muvec,rf)


    

In [307]:

    

fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_title('Portfolio Analysis')
ax.set_ylabel('Return')
ax.set_xlabel('Volatility')

x,y = zip(*frontier)
line, = ax.plot(x, y, color='b', lw=1)

ax.scatter(SD, Muvec, color='c')
ax.scatter(Equal_vol, Equal_return, color='k')

ax.scatter(0, rf, color='m')
ax.scatter((max(Muvec)-rf)/Sharpe_ratio, max(Muvec), color='m')
ax.scatter(max(SD), rf+max(SD)*Sharpe_ratio, color='m')

ax.plot([0,max(SD)], [rf, rf+max(SD)*Sharpe_ratio], color='y')

ax.scatter(XSh_SD, XSh_return, color='red')
ax.scatter(Xmin_SD, Xmin_return, color='green')
plt.show()


    

In [ ]: