使用Cplex求解均值方差模型

   在进行多目标进化算法求解投资组合问题时，由于多目标进化算法求得的是一组近似最优解，对于最大化收益和最小化风险的两目标M-V模型，需要找到一个基准，方便作为参考。这里利用Cplex工具求解均值-方差模型。

1、模型构建

   马科维茨的均值-方差模型，实质上是一个非线性的双目标的最优化问题，为的是达到最大化收益和最小化风险双目标的均衡状态。模型的数学描述如式所示：

   其中，其中V为投资组合的风险，希望其尽可能小。R为投资组合的收益，希望其尽可能大。X为投资到各个证券上的资产比例组成的权重向量。
双目标问题聚合为单目标问题：
   下面我们将该双目标问题转化为单目标优化问题，它是假设在投资组合收益率相同的情况下，把风险降为最小的模型。即仅保留最小化风险这个目标，将收益转化为约束。

2、数据准备

   数据准备对模型构建至关重要，原始数据集为100只股票的周盘价，通过统计计算，我们得到98只股票的期望收益率以及收益率的协方差矩阵。

3、Cplex求解

（1）首先导入所需函数库

import cplex
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

(2)数据准备

cov_matrix =pd.read_csv(r'data\股票收盘价集合\股票收益率矩阵.csv')
expect_return = pd.read_csv(r'data\股票收盘价集合\股票期望收益率.csv')
cov_matrix = cov_matrix.to_numpy()
expect_return = expect_return['0'].tolist()
return_max = max(expect_return)
return_min = min(expect_return)
StockNames = range(0,98)

(3)构建Cplex模型

def setproblemdata(p):
    p.objective.set_sense(p.objective.sense.minimize)    
    names = ['w' + str(i) for i in range(98)]    
    p.variables.add(obj=[0] * 98, ub=[1] * 98, lb=[0] * 98, names=names)
    my_senses = 'EG'
    my_rhs = [1, 0.005]
    rows = [[names, [1] * 98],[names, list(expect_return[:98])]]
    p.linear_constraints.add(lin_expr=rows, senses=my_senses, rhs=my_rhs )
    qmat = []
    for j in range(98):
        qmat.append([[i for i in range(98)], [2 * cov_matrix[j, m] for m in range(98)]])
    p.objective.set_quadratic(qmat)

（4）模型运行求解

p = cplex.Cplex()
setproblemdata(p)
p.solve()

print(sum(p.solution.get_values()))#是否满足预算约束
print(p.solution.get_objective_value())#最低风险值
print(np.array(p.solution.get_values()).round(3))#对应的投资组合

portfolio_names = []
portfolio_weights = []
for i in range(98):
    proportion = np.round(p.solution.get_values()[i], 3)  # keep results with 3 decimals and keep weights larger than 0
    if proportion != 0:
        portfolio_names.append(StockNames[i])
        portfolio_weights.append(proportion)
df = pd.DataFrame(columns = portfolio_names)
df.loc['Weight'] = portfolio_weights
print('Selected stocks and their weights:')
print(df)

运行结果

（5）敏感性分析

def sensitivity(exp_return):
    p = cplex.Cplex()
    p.objective.set_sense(p.objective.sense.minimize)
    names = ['w' + str(i) for i in range(98)]
    p.variables.add(obj=[0] * 98, ub=[1] * 98, lb=[0] * 98, names=names)
    my_senses = 'EG'
    my_rhs = [1, exp_return]
    rows = [[names, [1] * 98], [names, list(expect_return[:98])]]
    p.linear_constraints.add(lin_expr=rows, senses=my_senses, rhs=my_rhs)
    qmat = []
    for j in range(98):
        qmat.append([[i for i in range(98)], [2 * cov_matrix[j, m] for m in range(98)]])
    p.objective.set_quadratic(qmat)
    p.solve()
    return (p.solution.get_objective_value(), np.array(p.solution.get_values()).round(3))

return_values = np.arange(0.0000635, 0.0095, 0.00009436)
risks = []
for i in return_values:
    risks.append(sensitivity(i)[0])
print(return_values)
print(risks)
plt.figure(figsize=(15,8))
plt.plot(risks, return_values, '.--')
plt.xlabel('Portfolio Risk', fontsize = 18)
plt.ylabel('Portfolio Return', fontsize = 18)
plt.title('Efficient Frontier for our stocks pool\n', fontsize = 25)
plt.grid()
plt.show()
print(np.sum(sensitivity(0.005)[1]))

最优边界图片