使用Python描绘Markowitz有效边界
0 准备工作
在拙作《使用pandas-datareader包下载雅虎财经股价数据》中，我们使用pandas-datareader功能包对指定ticker和日期的股价数据进行下载，并整合了股利数据。本篇文章将使用pandas和numpy包对我们下载的数据进行处理，计算收益率，使用scipy包进行规划求解，并使用matplotlib包刻画Markowitz有效边界。
首先，我们要知道Markowitz有效边界到底是什么东西。从纯应用的角度看，这个有效边界是一条曲线，在坐标系中，横坐标是资产组合的风险（用收益率标准差σ表示），纵坐标是资产组合的收益率R。在这个坐标平面中，每一个点代表一个资产组合，它们拥有各自的收益率和风险水平。有效边界上的点，是针对同一收益率水平，可行域（指只使用给定universe中的资产，权重和小于等于100%的区域）中风险最小的资产组合。我们其实可以猜一猜这条有效边界的形状，大概就是风险越大收益越大的那种增函数的感觉。有一个大概的印象就好，因为我们就要把它画出来。
其次，我们要对计算收益率的假设有一定的认知。Yahoo Finance可以提供指定证券的adjusted closing price. 我们假设证券所发放的股利仍再投资于原证券，则Adjusted Closing Price的变化率可以反映该证券的真实收益率和真实波动率。在这里，我们使用对数收益率及其标准差。
对于给定universe中的资产， 我们还要计算它们的相关系数或协方差，这样，对于一个资产组合我们才能正确计算它的收益率和标准差。我们使用以下公式计算。
对于给定权重w，在收益率均值和标准差已知的情况下，我们可以计算某一投资组合的收益率均值和标准差。
最后，要掌握一点点线性代数的知识，因为收益率的计算十分容易理解，只要分别将权重和收益率相乘，再将各个乘积相加即可。而投资组合方差的两个求和公式展开以后要写很久。但学过线性代数的我们知道，组合方差可以用矩阵乘法求出。现在我们给定一个w.
那么，组合的方差可以由以下公式得出。
好神奇，我也不知道为什么（其实是线性代数），展开写果然就和上面的求和公式一模一样。
这样，我们就已经掌握了计算某一投资组合收益率和风险的方法，并可以根据收益率和风险，将其描绘在一个坐标平面上。完成这次作业，除了pands和pandas-datareader，我们还需要使用用于数值分析的功能包numpy，用于绘图的功能包matplotlib，以及可以用来做最优化（类似excel中solver功能）的功能包scipy，使用pip方法将其安装到你的Python中。现在，想想你最喜欢的股票，打开你的Python代码编辑器，你便完成了所有的准备工作。
1 最大限度获取有效历史数据
我们选择金融学教授们最喜欢的几只美股，Apple，Facebook，Google和Microsoft，作为我们的Universe，并将这些公司的名称储存在列表tickers中，公司的个数储存在变量q中。
import numpy as np
import pandas as pd
import pandas_datareader.data as web
import matplotlib.pyplot as plt
import scipy.optimize as solver
import datetime as dt
tickers = ['AAPL', 'FB', 'GOOG', 'MSFT']
q = len(tickers)
接下来，通过在《使用pandas-datareader包下载雅虎财经股价数据》中介绍过的方法，将下载到的2010年1月1日到今天的“Adjusted Closing Price”存储到一个新建的名为database的DataFrame中。这个部分的核心代码如下。
database = pd.DataFrame()
for i in tickers:
database[i] = web.DataReader(i, 'yahoo', '1/1/2010', datetime.date.today())['Adj Close']
考虑到Facebook上市的时间不过5年，在2010年1月1日时，它还没有股价，我们需要引入一个变量records，用它记录我们能使用的最大日期间隔，这个变量我们将会在后面计算收益率时用到。改良后的代码如下：
import pandas as pd
import numpy as np
import pandas_datareader.data as web
import matplotlib.pyplot as plt
import scipy.optimize as solver
import datetime as dt
database = pd.DataFrame()
tickers = ['AAPL', 'FB', 'GOOG', 'MSFT']
q = len(tickers)
records = 0
for i in tickers:
tmp = web.DataReader(i, 'yahoo', '1/1/2010', dt.date.today())
database[i] = tmp['Adj Close']
if records == 0:
records = len(tmp)
else:
records = min(records, len(tmp))
到此为止，我们已经下载到了我们需要的股价数据。
2 计算单个资产年化收益率和标准差
首先来计算单个资产的每一天的收益率。每计算一个日收益率，我们都要使用到“今天”和“昨天”两天的数据。在这里，使用DataFrame类型的一个Attribute，.shift(), 来实现收益率的计算。database.shift(1)表示database这个DataFrame所有数据往下平移一行，自然，其第一行变成了空值NaN.
要注意，Python本身并没有求对数的功能，因此要调用numpy包中的log函数，在之后的计算中，我们还会用到numpy包中的sqrt函数来开平方（当然也可以使用** .5来实现开平方）。
returns = np.log(database / database.shift(1))
returns这个DataFrame有很多缺陷，比如有些日子有些证券停牌了，没有发布股价，因此当天的return是NaN，再比如说Facebook在2010年还没有上市，那段时间的股价和收益率也是NaN，因此，我们要剪除未上市的记录，并将上市后停牌的日期的收益率填补为0， 为了实现这两个目的，我们用到了DataFrame的两个Attribute，一个是.tail(), 它的功能是返回DataFrame最后的指定条数的数据；另一个是.fillna()，它的作用是为了填补该DataFrame中的NaN，将其填补为我们需要的值，并返回新的DataFrame.
考虑到Facebook上市当天是没有收益率的（没有前一天的股价），因此我们只对records - 1条数据感兴趣。因此，这样进行修剪，并打印returns的前几行和后几行。
returns = returns.tail(records - 1)
returns.fillna(value=0, inplace=True)
print returns.head()
print returns.tail()
打印的结果如下。
AAPL FB GOOG MSFT
Date
2012-05-21 0.056626 -0.116378 0.022578 0.016266
2012-05-22 -0.007708 -0.093255 -0.021912 0.000336
2012-05-23 0.024107 0.031749 0.014311 -0.022083
2012-05-24 -0.009226 0.031680 -0.009562 -0.001375
2012-05-25 -0.005374 -0.034497 -0.020299 -0.000344
AAPL FB GOOG MSFT
Date
2017-05-08 0.026825 0.005443 0.007704 -0.000870
2017-05-09 0.006384 -0.003847 -0.002282 0.001449
2017-05-10 -0.004752 -0.001263 -0.003643 0.003903
2017-05-11 0.008611 -0.001665 0.001958 -0.012340
2017-05-12 0.013869 0.001931 0.001739 -0.001169
语句returns.fillna(value=0, inplace=True)中参数inplace我们给定了True，如果不这么做的话，该语句将只返回新的DataFrame，却并不改变旧的DataFrame. 用更易于理解的话说，该语句和下面语句的效果是一样的。
returns = returns.fillna(value=0)
利用DataFrame类型自带的.mean()和.cov()两个attribute，可以直接返回这四个证券的均值（一个pandas Series）和协方差矩阵（一个pandas DataFrame）。我们对返回后的值进行年化调整并打印。
mean = returns.mean() * 252
cov = returns.cov() * 252
print mean
print cov
打印结果如下。
AAPL 0.165849
FB 0.275372
GOOG 0.228091
MSFT 0.197116
dtype: float64
AAPL FB GOOG MSFT
AAPL 0.063315 0.020556 0.017735 0.019051
FB 0.020556 0.145224 0.028971 0.019541
GOOG 0.017735 0.028971 0.049977 0.024002
MSFT 0.019051 0.019541 0.024002 0.052169
Bravo! pandas功能包的descriptive stats功能非常简洁好用，直接给出了我们想要的均值序列和协方差矩阵，所有单个资产的收益率均值、方差（及标准差），以及相关性都包含在了mean和cov这两个变量中。
3 计算某一投资组合的年化收益率和风险水平
为了刻画有效边界，我们要假定市场不能卖空。如果出现了卖空，组合可能有很大的标准差，离原点很远，比较难以观测。
首先，我们随便写一个权重组合w，使之相加等于1. 比如0.1、0.2、0.3和0.4. 在给权重组合w赋值时，将其由列表转换为numpy功能包中的array类型，这样做是因为之后的计算中要用到w的转置矩阵（列向量），而list本身不能被转置。
w = np.array([.1, .2, .3, .4])
在进行下一步之前，要先普及两个函数，numpy.dot()和reduce().
在Python的矩阵运算中，*可以用来表示点乘，即对应位置乘法，我们可以用它来计算组合的收益率水平。而numpy包中的dot()函数的主要作用是用来做矩阵乘法，它的参数中主要的数据类型是numpy中的array类型。非常可喜的是，DataFrame类型也可以作为dot()函数的参数，也就是说，我们不用对cov变量做变成array类型的转换，即可让其参加dot()运算。
其实，numpy中还有一个matrix类型，两个matrix类型之间可以直接用*做矩阵乘法。我将在之后一篇介绍Value at Risk的笔记中讨论这个类型的用法。
我们用w.T表示w的转置，则计算组合收益率和标准差可以分成以下三步：
r = sum(w * mean)
var = np.dot(w, cov)
var = np.dot(sds, w.T)
s = np.sqrt(var) # s = var ** .5
如何简化这个步骤，而并不使第一行中的var产生歧义呢？dot()函数是个二元函数，我们可以使用Python内置的reduce()函数来简化这个过程。reduce()函数一般有两个参数，第一个参数是二元函数（也可以是支持二元情况的多元函数，如min函数）的名字，本例中为np.dot, 另外一个参数是一个list，它有两个或以上个元素。运行时，reduce函数会先用list中的前两个元素进行指定的函数运算，再用运算结果和第三个元素进行运算，依次类推。最后返回的结果，就是所有元素都参与了运算的最终结果。
有了这个函数，我们可以舍弃var这个变量，简化之前的计算方法。
r = sum(w * mean) # multiply
s = np.sqrt(reduce(np.dot, [w, cov, w.T])) # dot multiply
获取了第一个投资组合的收益率和标准差，我们已经迫不及待地想把它表示在坐标平面上看看它是什么样子了。为了将它描绘在坐标平面上，我们在一开始引用了matplotlib包中的pyplot模块。现在要调用该模块中的plot方法。最直观的两个参数，自然就是x和y了。为了纪念它，我们可以用一颗小星星把它标出来，将第三个参数改为'y*', 代表yellow star。
plt.plot(r, s, 'y*')
如果你是在Console或Jupyter Notebook中键入了以上命令，下面这张图应该会直接弹出或显示出来。
夜空中最亮的星
如果你使用PyCharm来编写代码，也很棒，在PyCharm的Tools菜单中，有一个Python Console命令，可以呼出Console，非常方便。同时，PyCharm也支持Jupyter Notebook，可以新建ipython文件。
Python Console
4* 使用蒙特卡洛方法找出可行域
这一部分是一个可选部分，不会影响最后的结果。
我们知道，可行域中有无数的投资组合，我们如果只用拍脑门的方法，寻找诸如[.1, .2, .3, .4]这样的权重组合，效率很低。应该找一种更加系统的方法。《Python for Finance: Analyze Big Financial Data》中介绍了一种较为系统的生成权重的方法，即蒙特卡洛方法。
简单来说，蒙特卡洛模拟就是根据某一函数不同参数的分布及其相关性，大量随机取样以模拟函数值。在这里，“某一函数”的“函数值”就是收益率和标准差，而“不同参数”就是q个不同的权重。
由于我们生成的权重之间相互独立，又不存在卖空问题，我们只需要用numpy.random模块中的rand()函数即可。rand()函数的用法，和excel中很像，生成[0, 1]闭区间中的随机数，取到每一个数的概率都相等（uniform分布）。rand()函数默认参数是1，即只返回一个随机数，现在我们用它返回q个随机数组成的list，使用rand(q)即可。
由于q个随机数相加并不等于1，我们将每个随机数除以list的总和，这样得到的就是一组有意义的权重。因为Markowitz有效边界是有效的，可以理解为手中头寸达到了full-employed的状态，所以我们只对可行域中权重加和为1的点感兴趣。
为了更好地描绘可行域，模拟的次数自然越大越好。我们选择100000次作为模拟次数，将每一次模拟的收益率和标准差，分别储存在rtn和sds两个列表中。我们使用plot()函数来看一看模拟的分布情况。
sds = []
rtn = []
for _ in range(100000):
w = np.random.rand(q)
w /= sum(w)
rtn.append(sum(mean * w))
sds.append(np.sqrt(reduce(np.dot, [w, cov, w.T])))
plt.plot(sds, rtn, 'ro') # ro for red dot
结果如下。
可行域，权重和为100%
观察这浩瀚星海，与你一开始想的不同，你发现原来还有那种明明风险很大但收益率却很低的辣鸡投资组合。如果使用的是Console，先不要关闭弹出的图像窗口。
5 描出有效边界
Markowitz有效边界其实已经很明显了，就是可行域的左边界，它代表了同一收益率水平下，风险水平最低的投资组合。但是，我们很难把蒙特卡洛中这些不连续的点给揪出来。因此，我们可以根据其定义，使用scipy.optimize模块中的minimize函数来最小化给定收益率水平下的标准差。
这个函数有许多必选或可选参数，首先最重要的是一个目标函数。这个函数的函数值就是我们最小化的目标。如果我们想实现最大化，只要最小化其负值即可。
根据这道题，我们为投资组合w量身定做一个函数，使其返回值为该投资组合的标准差。
def sd(w):
return np.sqrt(reduce(np.dot, [w, cov, w.T]))
该函数只读一个参数，即w，它的类型应是numpy中的array类型。
我们来分析一下这个sd这个函数。它的参数是一个array，它的限制条件有两个，首先，w的加和为1；其次，w都是非负数。
minimize函数的参数fun就是我们要优化的函数的名字sd. 参数x0是一个初始值，正如excel中solver运行之前，你要给参数单元格赋一个初始值一个道理，这里，我们不妨赋给各个证券等权重。
x0 = np.array([1.0 / q for x in range(q)])
该方法可以生成一组长度为q，所有元素都是1/q的数组，并赋值给x0.
下一个我们关心的参数是bounds，它代表了目标函数输入参数的取值范围，在这里我们的输入的参数是权重array，它要求每一个权重都大于0. bounds参数的类型是tuple，这个类型我将在以后的笔记中专门介绍。如果输入为一个变量，那么这里用一个tuple(min, max)来规定其上界和下界，如果输入为一个list或array，那么用多个(min, max)组成的tuple((min, max), (min, max), ...)来规定每一个数值的上界和下界。这里，每一个权重的取值范围都是0到1.
bounds = tuple((0, 1) for x in range(q))
最后一个我们比较关心的参数就是constraints了。我们输入的权重，它们有两个constraints. 首先，它们的加和为1；其次，它们的收益率是我们给定的水平。constraints参数的类型是一个字典组成的list，list中的每一个字典一般都有两个key，其中一个叫做'type', 另一个叫做'fun', type的值可以是'eq'或者'ineq', 而fun的值是一个关于输入值的函数，用lambda表达式来表示。如果type为eq，则约束fun的函数值等于0，如果type为ineq，则约束fun的函数值大于0. 假设我们给定的收益率水平是0.18.
constraints = [{'type': 'eq', 'fun': lambda x: sum(x) - 1},
{'type': 'eq', 'fun': lambda x: sum(x * mean) - .18}]
确定了所有的参数，我们来试着运行一下minimize函数，将其结果保存到变量outcome中。
def sd(w):
return np.sqrt(reduce(np.dot, [w, cov, w.T]))
x0 = np.array([1.0 / q for x in range(q)])
bounds = tuple((0, 1) for x in range(q))
constraints = [{'type': 'eq', 'fun': lambda x: sum(x) - 1},
{'type': 'eq', 'fun': lambda x: sum(x * mean) - .18}]
outcome = solver.minimize(sd, x0=x0, constraints=constraints, bounds=bounds)
print outcome
打印结果如下。
fun: 0.19747746355286486
jac: array([ 0.22423773, 0.10300811, 0.10652425, 0.16511219])
message: 'Optimization terminated successfully.'
nfev: 30
nit: 5
njev: 5
status: 0
success: True
x: array([ 0.57070503, 0. , 0.02351972, 0.40577525])
简单地读一下结果，我们发现我们成功了！我们比较感兴趣的是结果中的fun和x两个值，他们分别是目标函数的最小值（min(sd)）和对应的参数输入（w*）。我们可以把它们当成outcome的两个Attribute进行调用。不妨来验算一下这个投资组合的收益率是否是0.18.
print sum(outcome.x * mean)
结果如下。
0.18000000000191341
非常接近0.18. 接下来我们要修改刚才的优化准备过程，写一个循环语句，将每一个给定水平的最低风险储存到一个列表中。在这之前，先建立一个给定收益率array，使用numpy中的arange()函数，设定所有给定收益率在图中显示的0.18到0.26之间，每隔0.005记数一次。
given_r = np.arange(.18, .26, .005)
arange()函数的好处是初值、终值、步长皆可是小数，而python自带的range函数的三个参数只能是整数。
接下来，将我们的优化过程改写为循环语句。
def sd(w):
return np.sqrt(reduce(np.dot, [w, cov, w.T]))
x0 = np.array([1.0 / q for x in range(q)])
bounds = tuple((0, 1) for x in range(q))
given_r = np.arange(.16, .28, .005)
risk = []
for i in given_r:
constraints = [{'type': 'eq', 'fun': lambda x: sum(x) - 1},
{'type': 'eq', 'fun': lambda x: sum(x * mean) - i}]
outcome = solver.minimize(sd, x0=x0, constraints=constraints, bounds=bounds)
risk.append(outcome.fun)
现在，我们要在原来浩瀚星海的基础上描有效边界了。如果刚刚不小心把浩瀚星海给关了，只要重新做一遍plot命令即可。现在，将有效前沿描出来。
plt.plot(risk, given_r, 'x')
出来大概是这个样子的。
浩瀚星海的边界
如果我们想找出全局最小标准差所在的位置，只需将优化的constraints中的给定收益率一条删除即可，因为在这浩瀚星海（universe）中，只有一颗星星拥有最小的标准差，而它的收益率水平也是唯一的。
constraints = {'type': 'eq', 'fun': lambda x: sum(x) - 1}
minv = solver.minimize(sd, x0=x0, constraints=constraints, bounds=bounds).fun
minvr = sum(solver.minimize(sd, x0=x0, constraints=constraints, bounds=bounds).x * mean)
plt.plot(minv, minvr, 'w*') # w* for white star
对我们做好的图表做一些修饰。
plt.grid(True)
plt.title('Efficient Frontier: AAPL, FB, GOOG, MSFT')
plt.xlabel('portfolio volatility')
plt.ylabel('portfolio return')
最后的结果如下。
Efficient Frontier: AAPL, FB, GOOG, MSFT
这篇笔记是FIN 567课程的最后一次作业的笔记，它是一个典型的金融学作业和Python练习题目，包含了股价的下载、收益率和协方差的计算，补全NaN值，矩阵的点乘和叉乘，蒙特卡洛模拟，最优化问题，Python绘图等。尤其是其中的minimize方法，它几乎可以解决所有excel中solver可以解决的问题，说它是Python中的solver也不为过。
虽然这个作业用到了蒙特卡洛模拟，但由于我们所模拟的权重实际上是确定的值，而且模拟的数值之间并没有关系，因此可以说与蒙特卡洛模拟的关系并不大。今后将会有一篇笔记记录使用Python实现蒙特卡洛模拟其他例子。
关于Python在金融学中的应用，还是要推荐一本教材，《Python for Finance: Analyze Big Financial Data》。这本书成书较早，其中用到的一些函数如今已经有了更新，但是算法部分还是非常值得学习的。
此外，到现在为止，我还在使用Python 2，但逐渐想向Python 3靠拢，毕竟Python 3今后是要取代Python 2的。所以今后的笔记可能会依靠Python 3的语法。如果您发现任何问题或有任何疑问，欢迎指正或讨论。
by JohnnyMOON, COB @UIUC
这只是做Finance作业的学习笔记
EM: gengyug2@illinois.edu
