Solution

def portfolioVariance(weights,matrix):
    weights = np.mat(weights)
    matrix = np.mat(matrix)
    return (weights*matrix*weights.transpose()).item()*252
portfolioVariance([.6,.2,.2],df.cov())

We need to use item() to get the item out of the matrix!

I am going to also present a secondary way to do this, for your knowledge. The below function sums all the variances for each stock adjusted for weight, and then finds each pairwise covariance and adds it with a weight adjustment multiplied by 2. The reason we are going to use the first method, however, is because it simply is more effecient.

def portfolioVariance2(weights,matrix):
    pvar = 0
    for i in range(len(weights)):
        var = matrix.iloc[i,i]
        pvar+=weights[i]**2*var
    for i1 in range(len(weights)):
        for i2 in range(i1+1,len(weights)):
            var1 = matrix.iloc[i1,i1]
            var2 = matrix.iloc[i2,i2]
            cov = matrix.iloc[i1,i2]
            pvar+=2*weights[i1]*weights[i2]*cov
    return pvar*252
portfolioVariance2([.6,.2,.2],df.cov())

We are going to need a portfolio return function as well.

def portfolioReturn(weights,returns):
    return sum([x*y for x,y in zip(weights,returns)])

rets = [.09,.08,.05]
portfolioReturn([.2,.3,.5],rets)

Here comes the fun part, we are going to get every combination of the three weights for stocks (incrementing by .01), and see what the different opportunities look like for investing in different weights.

r = []
v = []
for x in range(101):
    for y in range(0,101-x):
        z=100-x-y
        r.append(portfolioReturn([x/100,y/100,z/100],rets))
        v.append(portfolioVariance([x/100,y/100,z/100],df.cov()))
std = [x**.5 for x in v]

The first two lines are out return and variance arrays being created. Next, we loop so that we get every possible combo. In our y loop, we end it at 101-x because we can’t have an x and y value that add up to more than 100%. Finally, z is whatever weight is leftover. The rest of the loop is running the functions to get the values.

import matplotlib.pyplot as plt
plt.plot(std,r)
plt.xlabel("Standard Deviation")
plt.ylabel("Expected Return")
plt.title("Portfolio Sets")
plt.show()

And here is what we can get from combining the stocks in different ways! Notice the effecient frontier is still the most outer points along the shape, where there are no points with the same risk and higher return (they are the highest return portfolio for that level of risk).

Add in some points representing 100% weight in one stock.

import matplotlib.pyplot as plt
plt.plot(std,r)
plt.xlabel("Standard Deviation")
plt.ylabel("Expected Return")
plt.title("Portfolio Sets")
for i in range(len(rets)):
    plt.plot((df.cov().iloc[i,i]*252)**.5,rets[i],"ro")
plt.show()

Let’s cover how to make an array of all ones first.

print(np.ones(10))

The argument defines the length of the array.

Now, let’s say we want 10 values that add up to 1, and we want them to be random. If we feed np.random.dirichlet() an array of all 1’s then we can get back a random set of decimals that sum to 1, we will see why this is important very soon.

import numpy as np
print(np.random.dirichlet(np.ones(10)))

Why would we want to do this? Well, we want to simulate what we did before but with 10 stocks. Trying to compute this many different possible weights would take a massive amount of time, so instead we are going to do a monte carlo simulation. What we are going to do is test a lot of different random sets of weights, and after we sample a large amount we expect that it will look somewhat like the real distribution of every possible combination.

Time to get the data for 10 stocks.

import pandas_datareader.data as web
import datetime
start = datetime.datetime(2016, 1, 1)
end = datetime.datetime(2017, 1, 1)
df = pdr.get_data_yahoo(stocks[10:20], start, end)["Adj Close"]
df = df.pct_change().dropna()

For expected return, we will use 2*standard deviation+.02. This isn’t a measure you would ever use in reality, but it helps to have expected returns that reflect risk for modeling this.

print((np.std(df)*2+.02))

We can see these returns look reasonable. Let’s make them into an array by using tolist().

rets = (np.std(df)*2+.02).tolist()
matrix = df.cov()