1. 问题描述
现假设有A, B, C, D, E五只股票的收益率数据((第二日收盘价-第一日收盘价)/第一日收盘价)), 如果投资人的目标是达到20%的年收益率,那么该如何进行资产配置,才能使得投资的风险最低?
更一般的问题,假设现有x1,x2,...,xn, n支风险资产,且收益率已知,如果投资人的预期收益为goalRet,那么该如何进行资产配置,才能使得投资的风险最低?
2. Markowitz 投资组合理论简介
1952年，芝加哥大学的Markowitz提出现代资产组合理论(Modern Portfolio Theory，简称MPT)，为现代西方证券投资理论奠定了基础。其基本思想是，证券投资的风险在于证券投资收益的不确定性。如果将收益率视为一个数学上的随机变量的话，证券的期望收益是该随机变量的数学期望（均值），而风险可以用该随机变量的方差来表示。
对于投资组合而言，如何分配各种证券上的投资比例，从而使风险最小而收益最大？
答案是将投资比例设定为变量，通过数学规划，对每一固定收益率求最小方差，对每一个固定的方差求最大收益率，这个多元方程的解可以决定一条曲线，这条曲线上的每一个点都对应着最优投资组合，即在给定风险水平下，收益率最大，这条曲线称作“有效前沿” (Efficient Frontier)。
对投资者而言，不存在比有效前沿更优的投资组合，只需要根据自己的风险偏好在有效前沿上寻找最优策略。
简化后的公式为:
其中 p 为投资人的投资目标,即投资人期待的投资组合的期望值. 目标函数说明投资人资产分配的原则是在达成投资目标
p的前提下,要将资产组合的风险最小化,这个公式就是Markowitz在1952年发表的'Portfolio Selection'一文的精髓,该文奠定了现代投资组合理论的基础,也为Markowitz赢得了1990年的诺贝尔经济学奖. 公式(1)中的决策变量为wi, i = 1,...,N, 整个数学形式是二次规划(Quadratic Programming)问题,在允许卖空的情况下(即wi可以为负,只有等式约束)时,可以用拉格朗日(Lagrange)方法求解。
有效前缘曲线如下图:
3. 二次规划求解方法介绍
3.1 等式约束凸二次规划的解法
我们考虑如下的二次规划问题
运用拉格朗日方法求解,可以得到
再看公式(1),则将目标函数由 min WTW 调整为 min 1/2(WT
W), 两问题等价,写出的求解矩阵为:
3.2 含有不等式约束的凸二次规划的解法
关于该问题的解法,可利用KKT条件来进行求解,具体的方法可以参考:
拉格朗日乘子法和KKT条件,这里给大家简单介绍一下python中求解二次规划问题的包 CVXOPT。
工具包: CVXOPT python凸优化包
函数原型: CVXOPT.solvers.qp(P,q,G,h,A,b)
求解时,将对应的P,q,G,h,A,b写出,带入求解函数即可.值得注意的是输入的矩阵必须使用CVXOPT 中的matrix函数转化,输出的结果要使用 print(CVXOPT.solvers.qp(P,q,G,h,A,b)['x']) 函数才能输出。
4. 实例
这里选取五支股票2014-01-01到2015-01-01的收益率数据进行分析.
选取的五支股票分别为: 白云机场, 华夏银行, 浙能电力, 福建高速, 生益科技
先大体了解一下五支股票的收益率情况:
看来，20%的预期收益是达不到了。
接下来,我们来看五支股票的相关系数矩阵：
可以看出，白云机场和福建高速的相关性较高，因为二者同属于交通版块。在资产配置时，不利于降低非系统性风险。
接下来编写一个ＭeanVariance类，对于传入的收益率数据，可以进行给定预期收益的最佳持仓配比求解以及有效前缘曲线的绘制。
绘制的有效前缘曲线为：
将数据分为训练集和测试集，并将随机模拟的资产配比求得的累计收益与测试集的数据进行对比，得到：
可以看出，在前半段大部分时间用Markowitz模型计算出的收益率要高于随机模拟的组合，然而在后半段却不如随机模拟的数据，可能是训练的数据不够或者没有动态调仓造成的，在后面写策略的时候，我会加入动态调仓的部分。
5. 代码
股票分析部分：
# -*- coding: utf-8 -*-
# @Time : 2019/2/9 0:05
# @Author : Arron Zhang
# @Email : 549144697@qq.com
# @File : plot return rate.py
# @Software: PyCharm
import pandas as pd
import baostock as bs
import matplotlib
from matplotlib import pyplot as plt
import numpy as np
#获取给定时间段的股票交易信息
def get_stock_data(t1,t2,stock_name):
lg = bs.login()
print('login respond error_code:' + lg.error_code)
print('login respond error_msg:' + lg.error_msg)
#### 获取沪深A股历史K线数据 ####
# 详细指标参数，参见“历史行情指标参数”章节
rs = bs.query_history_k_data(stock_name,
"date,code,open,high,low,close,preclose,volume,amount,adjustflag,turn,tradestatus,pctChg,isST",
start_date=t1, end_date=t2,
frequency="d", adjustflag="3")
print('query_history_k_data respond error_code:' + rs.error_code)
print('query_history_k_data respond error_msg:' + rs.error_msg)
#### 打印结果集 ####
data_list = []
while (rs.error_code == '0') & rs.next():
# 获取一条记录，将记录合并在一起
data_list.append(rs.get_row_data())
result = pd.DataFrame(data_list, columns=rs.fields)
print(result)
#### 结果集输出到csv文件 ####
result.to_csv("D:\stockdata\history_A_stock_k_data.csv", index=False)
print(result)
#### 登出系统 ####
bs.logout()
result['date'] = pd.to_datetime(result['date'])
result.set_index("date", inplace=True)
return result
byjc = get_stock_data('2014-1-1','2015-1-1','sh.600004')
hxyh = get_stock_data('2014-1-1','2015-1-1','sh.600015')
zndl = get_stock_data('2014-1-1','2015-1-1','sh.600023')
fjgs = get_stock_data('2014-1-1','2015-1-1','sh.600033')
sykj = get_stock_data('2014-1-1','2015-1-1','sh.600183')
by = byjc['pctChg']
by.name = 'byjc'
by = pd.DataFrame(by,dtype=np.float)/100
hx = hxyh['pctChg']
hx.name = 'hxyh'
hx = pd.DataFrame(hx,dtype=np.float)/100
zn = zndl['pctChg']
zn.name = 'zndl'
zn = pd.DataFrame(zn,dtype=np.float)/100
fj = fjgs['pctChg']
fj.name = 'fjgs'
fj = pd.DataFrame(fj,dtype=np.float)/100
sy = sykj['pctChg']
sy.name = 'sykj'
sy = pd.DataFrame(sy,dtype=np.float)/100
sh_return = pd.concat([by,fj,hx,sy,zn],axis=1)
#了解各股票的收益率状况,并作图
sh_return = sh_return.dropna()
cumreturn = (1+ sh_return).cumprod()
sh_return.plot()
plt.title('Daily Return of 5 Stocks(2014-2015)')
plt.legend(bbox_to_anchor = (0.5,-0.3),ncol = 5,fancybox = True,shadow = True)
cumreturn.plot()
plt.title('Cumulative Return of 5 stocks(2014-2015)')
Markowitz 投资组合模型求解
# -*- coding: utf-8 -*-
# @Time : 2019/2/11 11:10
# @Author : Arron Zhang
# @Email : 549144697@qq.com
# @File : Markowitz portfolio selection.py
# @Software: PyCharm
#构建一个MeanVariance类,该类可以根据输入的收益率序列,求解二次规划问题,计算出最优资产比例,并绘制最小方差前缘曲线
#定义 MeanVariance类
from matplotlib import pyplot as plt
from scipy import linalg
import numpy as np
import ffn
import cvxopt
from cvxopt import matrix
class MeanVariance:
#定义构造器,传入收益率数据(dataframe格式的每日收益率)
def __init__(self,returns):
self.returns = returns
#定义最小化方差的函数,即求解二次规划
def minVar(self,goalRet):
covs = np.array(self.returns.cov())
means = np.array(self.returns.mean())
L1 = np.append(np.append(covs.swapaxes(0,1),[means],axis=0),
[np.ones(len(means))],axis=0).swapaxes(0,1)
L2 = list(np.ones(len(means)))
L2.extend([0,0])
L3 = list(means)
L3.extend([0,0])
L4 = np.array([L2,L3])
L = np.append(L1,L4,axis=0)
results = linalg.solve(L,np.append(np.zeros(len(means)),[1,goalRet]))
return np.array([list(self.returns.columns), results[:-2]])
#定义绘制最小方差前缘曲线函数
def frontierCurve(self):
goals = [x/500000 for x in range(-100,4000)]
variances = list(map(lambda x: self.calVar(self.minVar(x)[1,:].astype(np.float)),goals))
plt.plot(variances,goals)
#给定各资产的比例,计算收益率的均值
def meanRet(self,fracs):
meanRisky = ffn.to_returns(self.returns).mean()
#不符合条件时,弹出错误
assert len(meanRisky) == len(fracs), 'length of fractions must be equal to number of assets'
return np.sum(np.multiply(meanRisky,np.array(fracs)))
#给定各资产的比例,计算收益率方差
def calVar(self,fracs):
#np.dot 可以将dataframe类型与矩阵直接相乘,得到的结果是array
return (np.dot(np.dot(fracs,self.returns.cov()),fracs))
#sim_weight = np.apply_along_axis(lambda x: x/sum(x),1,xim_weight)
#apply_along_axis可以将函数作用于矩阵的行或者列
#lambda函数的作用是对传入的参数直接给出结果,如果需要传入多个参数,可以写成array的形式传入,此案列中应用到了对数据的标准化中
def solve_quadratic_problem(self,goal):
covs = np.array(self.returns.cov())
means = np.array(self.returns.mean())
P = matrix(np.dot(2,covs))
Q = matrix(np.zeros((len(means),1)))
G = -matrix(np.zeros((len(means),len(means))))
A = matrix(np.append([np.ones(len(means))],[means],axis=0))
h = matrix(np.zeros((len(means),1)))
b = matrix(np.array([[1,goal]]).swapaxes(0,1))
sol = cvxopt.solvers.qp(P, Q, G , h, A, b)
return sol
6. 参考资料
蔡立耑：量化投资——以python为工具. 电子工业出版社
