``````

In [1]:

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from datetime import datetime

``````
``````

In [2]:

%matplotlib inline

``````
Run the next three cells for real data
``````

In [176]:

start = datetime(2004,2,1)
end = datetime.today()

# stock list
L = ['DDD', 'DIS', 'COST', 'SBUX']

#set up DataFrames
daily_price  = pd.DataFrame(index=pd.bdate_range(start, end)) # business days
daily_return = pd.DataFrame(index=pd.bdate_range(start, end))

``````
``````

In [177]:

# get daily equity "Adj Close" from start to end
# would like to build a database of SP500 stocks instead

daily_return = np.log(1+daily_price.pct_change())  # for a continuous return number
# cumulative_return = daily_return.cumsum()        useful for a normalized return chart

``````
``````

In [178]:

# create expected return vector, stdev vector and covariance, correlation matrices

R = daily_return.mean() # expected return vector
AAR = (1+R)**252-1      # average annual return vector
SD = daily_return.std()  # expected standard deviation vector
C = np.matrix(daily_return.cov())  # covariance matrix
Corr =  daily_return.corr() # and a correlation matrix for info

``````
Comment out R, SD and C in the cell below to run with real data from above
``````

In [215]:

W = [1/3, 1/3, 1/3, 0]
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 = 0.01  # risk free rate

``````
``````

In [216]:

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(R):         # W is a list [] of weights as real numbers
Xvec = np.ones(len(R))/len(R)
return Xvec

``````
``````

In [217]:

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

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

``````
``````

In [219]:

Muvec, Onevec, Xvec, LMtarget = Initialize(R,r)

``````
``````

In [220]:

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.item()

``````
``````

Out[220]:

(0.039999999999999994, 0.04472135954999579)

``````
``````

In [221]:

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

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

``````
``````

In [223]:

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)
return Xvec_target, target_SD

``````
``````

In [224]:

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

``````
``````

Out[224]:

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

``````
``````

In [225]:

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

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

``````
``````

Out[226]:

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

``````
``````

In [227]:

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

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

``````
``````

Out[228]:

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

``````
``````

In [229]:

def Find_Frontier(C, Muvec, rf):
frontier = []
frontier_weights = []
r_points = np.arange(min(Muvec)/1.2, max(Muvec)*1.2, (max(Muvec)-min(Muvec))/100)
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((r, Target_weights))
return frontier, frontier_weights

``````
``````

In [230]:

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

``````
``````

In [235]:

fig = plt.figure()
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='orange')
# ax.scatter((max(Muvec)-rf)/Sharpe_ratio, max(Muvec), color='orange')
ax.scatter(max(SD), rf+max(SD)*Sharpe_ratio, color='orange')

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

ax.scatter(XSh_SD, XSh_return, color='g')
ax.scatter(Xmin_SD, Xmin_return, color='r')

plt.grid(True)
plt.show()

``````
``````

``````
``````

In [233]:

fig = plt.figure()
ax.set_title('Portfolio Analysis')
ax.set_ylabel('weight')
ax.set_xlabel('Return')

r, W = zip(*frontier_weights)
y1=[]
y2=[]
y3=[]
y4=[]
for point in W:
y1.append(point[0].item())
y2.append(point[1].item())
y3.append(point[2].item())
#    y4.append(point[3].item())

plt.plot(r,y1, color='b')
plt.plot(r,y2, color='c')
plt.plot(r,y3, color='r')
#plt.plot(r,y4, color='m')

plt.plot([Xmin_return.item(), Xmin_return.item()], [0, 1], color='orange', lw=1)
plt.plot([XSh_return.item(), XSh_return.item()], [0, 1], color='g', lw=1)

# plt.grid(True)

``````
``````

Out[233]:

[<matplotlib.lines.Line2D at 0x10bc1f110>]

``````
``````

In [ ]:

``````