之前写过一篇。那个时候不是很懂……
Review：数量金融学（3）：Markowitz均值-方差模型
现在重新整理一下。
第一部分，生成各个资产组合的风险-收益特征。
先上效果图。
再上代码。
# -*- coding: utf-8 -*-import numpy as npimport matplotlib.pyplot as pltimport pandas as pdimport seaborn as sns# Random seednp.random.seed(123)## NUMBER OF ASSETSn_assets = 4## NUMBER OF OBSERVATIONSn_obs = 1000return_vec = np.random.randn(n_assets, n_obs) # print(return_vec.shape) # n_assets个标的（行），n_obs次观测值（列）plt.plot(return_vec.T, alpha=.4);plt.xlabel('time')plt.ylabel('returns')plt.savefig('random_ret.jpg')# 生成随机权重def rand_weights(n): ''' Produces n random weights that sum to 1 ''' k = np.random.rand(n) return k / sum(k)# 生成随机组合def random_portfolio(returns): ''' Returns the mean and standard deviation of returns for a random portfolio ''' # 每日收益率序列是根据“已知（或估计）的正态ret随机变量生成的” p = np.asmatrix(np.mean(returns, axis=1)) # 随机生成各资产权重 w = np.asmatrix(rand_weights(returns.shape[0])) C = np.asmatrix(np.cov(returns)) mu = w * p.T sigma = np.sqrt(w * C * w.T) # This recursion reduces outliers to keep plots pretty # 后期要作图，希望让 风险-收益 的点，出现在一个适合显示的范围 # （因为不限制的话，这些点可能会飞到很远的地方去） if sigma > 1: # （如果std>1，后面就不显示这个点） return random_portfolio(returns) # 重新计算 return mu, sigma# In[23]:# 这里仿真n_portfolios个资产组合n_portfolios = 1000means, stds = np.column_stack([ random_portfolio(return_vec) for _ in range(n_portfolios)])# In[25]:plt.plot(stds, means, 'o', markersize=5)plt.xlabel('std')plt.ylabel('mean')plt.title('Mean and standard deviation of returns of randomly generated portfolios')plt.savefig('risk_ret.jpg')
关于column_stack函数的说明。
# np.column_stack
# Signature: np.column_stack(tup)
# Docstring:
# Stack 1-D arrays as columns into a 2-D array.
# Take a sequence of 1-D arrays and stack them as columns
# to make a single 2-D array. 2-D arrays are stacked as-is,
# just like with `hstack`. 1-D arrays are turned into 2-D columns
# first.
# Parameters
# ----------
# tup : sequence of 1-D or 2-D arrays.
# Arrays to stack. All of them must have the same first dimension.
# Returns
# -------
# stacked : 2-D array
# The array formed by stacking the given arrays.
# See Also
# --------
# hstack, vstack, concatenate
# Examples
# --------
# >>> a = np.array((1,2,3))
# >>> b = np.array((2,3,4))
# >>> np.column_stack((a,b))
# array([[1, 2],
# [2, 3],
# [3, 4]])
与np.cov相关的一个函数，correof
设index_ret_df是收益率时间序列矩阵（n_asset种资产，每种资产各一列，第t行是这些资产第t期的收益率）
corroef_df = pd.DataFrame(np.corrcoef(index_ret_df,rowvar=0)) # rowvar的意思是，列与列之间求相关系数。
第二部分，优化资产组合。
这时可以用到优化工具箱。（scipy也可以，cvxopt也可以，但是我知道为什么的py3装cxvopt好费劲……好像要先装一堆什么东西，不是靠简单的pip install就能装上的）
定义优化函数。在既定的预期收益率下，求解使得资产组合收益率标准差（风险）最小的各资产权重。
有一些约束条件。要把方程整理为cvxopt的标准形式。（scipy也一样，一样需要输入一些标准化的形式）
# -*- coding: utf-8 -*-import numpy as npimport matplotlib.pyplot as pltimport cvxopt as optfrom cvxopt import blas, solversimport pandas as pdimport seaborn as snsTurn off progress printing solvers.options['show_progress'] = Falsedef optimal_portfolio(returns): n = len(returns) returns = np.asmatrix(returns) N = 100 mus = [10**(5.0 * t/N - 1.0) for t in range(N)] # 注意！这个mu不是我们要求的最小预期收益率mu！后面我会附上一些材料来解释。 # Convert to cvxopt matrices S = opt.matrix(np.cov(returns)) pbar = opt.matrix(np.mean(returns, axis=1)) # # Create constraint matrices G = -opt.matrix(np.eye(n)) # negative n x n identity matrix h = opt.matrix(0.0, (n ,1)) A = opt.matrix(1.0, (1, n)) b = opt.matrix(1.0) # Calculate efficient frontier weights using quadratic programming portfolios = [solvers.qp(mu*S, -pbar, G, h, A, b)['x'] for mu in mus] # 这个最优化的做法和我们的直觉有点不一样（后面我也会说） ## CALCULATE RISKS AND RETURNS FOR FRONTIER returns = [blas.dot(pbar, x) for x in portfolios] risks = [np.sqrt(blas.dot(x, S*x)) for x in portfolios] ## CALCULATE THE 2ND DEGREE POLYNOMIAL OF THE FRONTIER CURVE m1 = np.polyfit(returns, risks, 2) x1 = np.sqrt(m1[2] / m1[0]) # CALCULATE THE OPTIMAL PORTFOLIO wt = solvers.qp(opt.matrix(x1 * S), -pbar, G, h, A, b)['x'] return np.asarray(wt), returns, risks
调用求解最优组合的函数。
optweights, optreturns, optrisks = optimal_portfolio(return_vec)plt.plot(stds, means, 'o')plt.ylabel('mean')plt.xlabel('std')plt.plot(optrisks, optreturns, 'y-o')plt.plot(statistics(opts['x'].reshape(4,1))[1], statistics(opts['x'].reshape(4,1))[0], 'r^')
关于用cxvopt来求解Markowitz最优化问题的一些补充材料。
0. cvxopt最优化二次型的标准形式
1. 问题描述
2. 问题求解
直观一点，有效前沿可以这么求：
给定一个预期的收益，求在满足各种约束的条件下的各资产权重，让组合的标准差最小。
但根据上面的补充材料，可以把目标函数写为
f = -p'x + u x'Sx
这相当于是把向量目标函数标量化了。（1. 要最大化收益。 2. 要最小化风险）
其中u是一个可变参数。
当u=0的时候，最小化f可以代表最大化收益；当u=无穷大的时候，最小化f可以代表最小化风险。
所以，看起来这么诡异！
至于代码中给出的最优资产配置权重的计算，我暂时还不知道怎么解释！！
切点模型的话，好像据说和最大化Shapre Ratio等价，即(p'x - rf) / sqrt(x'Sx)。
但我不知道和代码里给出的
## CALCULATE THE 2ND DEGREE POLYNOMIAL OF THE FRONTIER CURVE m1 = np.polyfit(returns, risks, 2) x1 = np.sqrt(m1[2] / m1[0]) # CALCULATE THE OPTIMAL PORTFOLIO wt = solvers.qp(opt.matrix(x1 * S), -pbar, G, h, A, b)['x']是不是一个意思……
Refs.
Convex Optimization – Boyd and Vandenberghe.