资产组合理论与实战（3) - 性能优化

使用向量化以优化性能

上一篇文章中，我们使用蒙特卡洛方法，运行5000次，大约花了5.2秒。我们可以使用向量化运算方法来提升这个速度：

我们可以使用下面的算法，把部分计算移出了循环：

from empyrical import sharpe_ratio
import time
start = time.time()
num_ports = 5000
all_weights = np.random.rand(num_ports, len(stocks))

shape = (len(all_weights), -1)
all_weights = all_weights/all_weights.sum(axis=1).reshape(shape)

# got 241 daily summed returns over 5000 iters
all_port_returns = np.dot(all_weights,returns.T)

all_cov = np.cov(all_port_returns).diagonal()

# 计算波动率
all_vol = []
all_sharpe = []

annual_return = np.prod((1 + all_port_returns), axis=1) - 1

for i in range(num_ports):
    all_vol.append(np.sqrt(np.dot(all_weights[i].T, np.dot(all_cov[i], all_weights[i]))))
    all_sharpe.append(sharpe_ratio(all_port_returns[i]))

end = time.time()
print("计算用时", end-start)

# 绘图
all_vol = np.array(all_vol)
pos = np.argsort(all_vol)

plt.scatter(all_vol, annual_return, c=all_sharpe, cmap='RdYlBu')
plt.colorbar(label='Sharpe Ratio')

max_sharpe_pos = np.argmax(all_sharpe)
plt.scatter(all_vol[max_sharpe_pos], annual_return[max_sharpe_pos], c='red',s=80)

不过，在按行计算波动率和夏普率的时候，我不得不做出妥协，改用循环。在这种情况下（投资组合中只有4只标的），我们现在的执行时间需要0.65秒。提升了大概8倍的速度（初始为5.2s左右）。

PyPortfolioOpt的作者曾在一个示例中用到了这样的方法：

n_samples = 10000
w = np.random.dirichlet(np.ones(ef.n_assets), n_samples)
rets = w.dot(ef.expected_returns)
stds = np.sqrt(np.diag(w @ ef.cov_matrix @ w.T))
sharpes = rets / stds

计算速度非常之快（小于0.1秒），但需要对其中的ef.cov_matrix进行替换，感兴趣的同学可以研究下它的源码。

对于一个仅有4支标的的组合来说，这个速度（不到1秒）显然可以令人满意了。但我们也在上一篇文章中，提出一个问题，随着标的数的增加，使用蒙特卡洛模拟方法，我们需要暴力搜索的空间也随之变大。那么，究竟会增大多少呢，此时还可不可行？

我们先来看只有两支标的情况。假设资产权重的间隔是1%，即某个资产要么分配1%，要么分配2%，而不会分配1.05%这样的非整数权重。这样我们就可以算得，如果要把这些权重分布全部模拟到，搜索空间将是：

A: 1-100 ~ 101个选择项
B: 1个选择项

当A的权重选定以后，在权重和为1的约束下，B就只有一个选择。因此，总共是101次搜索。

当如果标的数量上升，我们为了模拟尽可能多的情况，一定会扩大模拟次数，因此总的运行时长迅速上涨。同时随着标的数的增加，我们要搜索的次数需要呈指数级增加，而不是按标的数线性增加。

因此，搜索次数是一个组合问题。假设存在四个标的，并不是简单的4×100就可以模拟出所有组合，而是存在：
$C_{99}^3$种分配方案。

对任意n个标的，每个标的的权重间隔按1%计算的话，这个搜索空间将是：

$C_{101-(n-1)}^{n-1}$次。如果我们有50支标的，则需要搜索最少1.998e+21亿次！

!!! tip
可以使用 math.comb来计算组合数。

所以，使用蒙特卡洛方法，会存在无法穷尽所有样本空间的可能性。我们必须采用数学的方法，来优化这个求解过程。

优化算法–scipy.optimize

寻找给定收益率下的最小波动率，或者给定波动率下的最大夏普率，这实际上是一类常见的优化问题，即：

我们可以使用 Scipy.optimize工具库来求解类似问题。

scipy.optimize包的主要功能

• SciPy Optimize提供最小化(或最大化)目标函数的功能，可能存在约束。
• 它包括解决非线性问题(支持局部和全局优化算法) ，线性规划，约束和非线性最小二乘，根寻找和曲线拟合。
• 当我们需要优化输入函数的参数时，scipy.optimize包含很多有用的方法，并可以处理不同种类的函数

它提供以下方法：

• minimize_scalar( )：使一个变量最小化
• minimize( )：使多个变量最小化
• curve_fit( )：寻找最优曲线，适合a set of data
• root_scalar( )：求一个变量的零点
• root( )：求多个变量的零点
• linprog( )：在线性等式和不等式的约束下，最小化线性目标函数

在这些方法中，我们主要介绍minimize_scalar和minimize。

minimize_scalar

有一类问题，比如解方程，它的解是一个数字，此时就可以用 minimize_scalar 方法来求解。

例如，我们的函数为：

$y=3x^4-2x+1$

minimize_scalar( )就是帮助我们找到函数最小值(y最小)时精确的坐标。通俗地说，就是解方程。

from scipy.optimize import minimize_scalar

def objective_function(x):
    return 3 * x ** 4 - 2 * x + 1
# 定义需要求解的函数

res = minimize_scalar(objective_function)
# 求出函数的最小值

print(res)

我们重点解释一下输出结果，这将在后面继续使用：

 message: 
          Optimization terminated successfully;
          The returned value satisfies the termination criteria
          (using xtol = 1.48e-08 )
 success: True
     fun: 0.17451818777634331
       x: 0.5503212087491959
     nit: 12
    nfev: 15

这里的fun是函数值，x是变量值，如果我们把x的值代入到公式 $3\ast x^4-2\ast x+1$，我们会得得到0.17,即fun的值。

要注意，不是所有情况都能进行优化求解的。比如，$y=x^3$没有最小值，如果我们对它求解，会导致overflowerror.

另外的情况是，存在多个最小值的情况，此时minimize_scalar不能保证找到函数的全局最小值。

我们以$y=x^4-x^2$为例，它存在多个最小值，minimize_scalar ()不能保证找到函数的全局最小值。我们尝试一下：

from scipy.optimize import minimize_scalar
# help(minimize)

def objective_function(x):
    return x**4 - x**2

res = minimize_scalar(objective_function)
res

我们将得到如下结果：

 message: 
          Optimization terminated successfully;
          The returned value satisfies the termination criteria
          (using xtol = 1.48e-08 )
 success: True
     fun: -0.24999999999999994
       x: 0.7071067853059209
     nit: 11
    nfev: 14

这样求解出来的方程的根是0.707。但我们知道，该方程至少还有一个-0.707的根，因此， minimize_scalar默认情况下，是无法找出全部的根的。但我们可能通过设定参数，来控制用于优化的求解器。一共有三种设定方式：

• brent:布伦特算法，也是该函数的默认算法
• golden:黄金分割算法，根据已有研究，该方式效果略差于brent
• bounded:布伦特算法的“有界实现”

brent和golden属于bracket：这是一个由两个或三个元素组成的序列，它提供了对最小区域边界的初始猜测。但是，这些求解器并不保证所找到的最小值将在此范围内。

bounded属于bounds：这是一个两个元素的序列，严格限制搜索区域的最小值。因此当最小值在已知范围内时，限制搜索区域才是有用的。这个概念我们在后面还将进一步接触。我们来看使用bounds方式的例子：

# 当我们采用bounds方式限定时
res = minimize_scalar(objective_function, method='bounded', bounds=(-1, 0))
res

现在，我们求解的x将落在区间(-1,0)之间，我们因此得到：

 message: Solution found.
 success: True
  status: 0
     fun: -0.24999999999998732
       x: -0.707106701474177
     nit: 10
    nfev: 10

这次我们得到了期望中的负数解。

不过，最佳资产组合的求解要比这个复杂，我们必须使用另一个方法，即minimize方法，具体优化算法则是最小顺序二乘法。

【系列文章目录】
资产组合理论及实战(1)-基本概念
资产组合理论及实战(2) - 蒙特卡洛方法