From discussions and code at https://class.coursera.org/fe-001/forum/thread?thread_id=476

meanvariance.py

#
# https://class.coursera.org/fe-001/forum/thread?thread_id=476
#
# https://www.dropbox.com/s/ghfzu5sthgb2gba/meanvariance.py
#
import numpy as np
from scipy import optimize
import unittest

class TestMeanVarianceMarket(unittest.TestCase):
    def setUp(self):
        self.cov_m = np.mat([[0.0010, 0.0013, -0.0006, -0.0007, 0.0001,  0.0001,  -0.0004, -0.0004],
                        [0.0013, 0.0073, -0.0013, -0.0006, -0.0022, -0.0010, 0.0014,  -0.0015],
                        [-0.0006, -0.0013, 0.0599, 0.0276, 0.0635, 0.0230, 0.0330, 0.0480],
                        [-0.0007, -0.0006, 0.0276, 0.0296, 0.0266, 0.0215, 0.0207, 0.0299],
                        [0.0001, -0.0022, 0.0635, 0.0266, 0.1025, 0.0427, 0.0399, 0.0660],
                        [0.0001, -0.0010, 0.0230, 0.0215, 0.0427, 0.0321, 0.0199, 0.0322],
                        [-0.0004, 0.0014, 0.0330, 0.0207, 0.0399, 0.0199, 0.0284, 0.0351],
                        [-0.0004, -0.0015, 0.0480, 0.0299, 0.0660, 0.0322, 0.0351, 0.0800]])
        self.mean_ret = np.mat([3.15, 1.75, -6.39, -2.86, -6.75, -0.54, -6.75, -5.26]).T
        self.interest = 1.5

    def test_return_of_portfolio(self):
        mv = MeanVarianceMarket(self.interest, self.mean_ret, self.cov_m)
        p = np.mat([5.6762, -0.5214, 2.5467, -1.7475, -1.6448, 4.2951, -8.6474, 1.0431])
        self.assertAlmostEqual(mv.return_of_portfolio(p), 67.358, 3)

    def test_volatility_of_portfolio(self):
        mv = MeanVarianceMarket(self.interest, self.mean_ret, self.cov_m)
        p = np.mat([5.6762, -0.5214, 2.5467, -1.7475, -1.6448, 4.2951, -8.6474, 1.0431])
        self.assertAlmostEqual(mv.volatility_of_portfolio(p), 99.657, 2)

    def test_sharpe_ratio_of_portfolio(self):
        mv = MeanVarianceMarket(self.interest, self.mean_ret, self.cov_m)
        p = np.mat([5.6762, -0.5214, 2.5467, -1.7475, -1.6448, 4.2951, -8.6474, 1.0431])
        self.assertAlmostEqual(mv.sharpe_ratio_of_portfolio(p), 0.661, 3)

    def test_volatility_of_minimum_variance_portfolio(self):
        mv = MeanVarianceMarket(self.interest, self.mean_ret, self.cov_m)
        m_p = mv.compute_minimum_variance_portfolio()
        self.assertAlmostEqual(mv.volatility_of_portfolio(m_p), 2.923, 3)

    def test_mean_return_sharpe_portfolio(self):
        mv = MeanVarianceMarket(self.interest, self.mean_ret, self.cov_m)
        s_p = mv.compute_sharpe_portfolio()
        self.assertAlmostEqual(mv.return_of_portfolio(s_p), 5.720, 3)

    def test_slope_capital_market_line(self):
        mv = MeanVarianceMarket(self.interest, self.mean_ret, self.cov_m)
        self.assertAlmostEqual(mv.slope_of_cml(), 0.766, 3)

class MeanVarianceMarket(object):
    """MeanVarianceMarket is implementation of Mean Variance Market model.
        interest_rate = float describing risk free rate
        mean_return is column vector (np.mat) of mean returns for assets
        cov_matrix is covarince matrix (np.mat) for assets
    """
    def __init__(self, interest_rate,  mean_return, cov_matrix):
        super(MeanVarianceMarket, self).__init__()
        self.interest_rate = interest_rate
        self.mean_return = mean_return
        self.cov_matrix = cov_matrix

    def return_of_portfolio(self, portfolio):
        """Portfolio is row vector (np.mat)"""
        #without risk free asset
        if len(portfolio.T) == len(self.mean_return):
            return float(portfolio * self.mean_return)
        #with risk free asset
        else:
            return float((portfolio[:,:-1] * self.mean_return) + (portfolio[:,-1] * self.interest_rate))

    def volatility_of_portfolio(self, portfolio):
        if len(portfolio.T) == len(self.mean_return):
            return float(100 * np.sqrt((portfolio * self.cov_matrix * portfolio.T)[0]))
        else:
            return float(100 * np.sqrt((portfolio[:,:-1] * self.cov_matrix * portfolio[:,:-1].T)[0]))

    def sharpe_ratio_of_portfolio(self, portfolio, benchmark=None):
        if benchmark==None:
            benchmark=self.interest_rate
        return (self.return_of_portfolio(portfolio)-benchmark)/self.volatility_of_portfolio(portfolio)

    def compute_minimum_variance_portfolio(self):
        initial_vect = np.array([1./len(self.mean_return) for _ in self.mean_return])
        def weight_constraint(portfolio):
            return np.array([np.sum(portfolio) - 1])
        def volatility(portfolio):
            return self.volatility_of_portfolio(np.mat(portfolio))

        cons = ({'type': 'eq', 'fun': weight_constraint})
        min_var_portfolio = optimize.minimize(volatility, initial_vect, constraints = cons, method = 'SLSQP') #options={'disp': True, 'iter':500, 'acc': 0.00001, 'epsilon': 0.000000001}
        return np.mat(min_var_portfolio.x)

    def compute_sharpe_portfolio(self):
        excess_returns = self.mean_return - self.interest_rate
        positions = excess_returns.T * np.linalg.inv(self.cov_matrix)
        sharpe_portfolio = positions / np.sum(positions)
        return sharpe_portfolio

    def slope_of_cml(self):
        return self.sharpe_ratio_of_portfolio(self.compute_sharpe_portfolio())

    def capm_return(self, beta):
        s_p = self.compute_sharpe_portfolio()
        mean_excess_return = self.return_of_portfolio(s_p) - mv_model.interest_rate
        asset_return = (mean_excess_return * beta) + self.interest_rate
        return asset_return

if __name__ == '__main__':
    unittest.main()