专栏首页量化投资与机器学习如何更稳健的计算组合最优权重(附代码)

如何更稳健的计算组合最优权重(附代码)

本期遴选论文 来源:SSRN 标题:A ROBUST ESTIMATOR OF THE EFFICIENT FRONTIER October 15, 2019 作者:Marcos López de Prado

今天分享的论文是Marcos López de Prado 2019年的论文《A ROBUST ESTIMATOR OF THE EFFICIENT FRONTIER》本文主要有两个创新点。

首先,提出了一种新的解决方法 ,称为嵌套聚类优化(NCO),该方法解决凸优化问题中噪声及复杂的信号结构引起的不稳定性。其次,作者还采用了蒙特卡罗模拟方法(Monte Carlo Optimization Selection, 以下简称为MCOS)对多种最优化算法产生的误差进行了评估(包括NCO),这样就可以根据评估的结果选择最稳健的优化模型。

不是一般性,假设有个系统有N个随机变量,他们的期望用向量

\mu

表示,协方差矩阵用

V

表示。目标是找到一个权重向量

w

使得系统的方差最小,即:

\begin{array}{c} \min _{\omega} \frac{1}{2} \omega^{\prime} V \omega \\ \text { s.t.: } \omega^{\prime} a=1 \end{array}

在金融领域,这就是一个典型的组合优化问题,当a为向量1是最优组合就是minimum variance portfolio。而当a为向量u时,最优组合就是夏普最大的组合。其解析解为:

\omega^{*}=\frac{V^{-1} a}{a^{\prime} V^{-1} a}

这类问题称为凸优化(CVO),为了简单起见,后面的所有讨论都基于这个最基本的凸优化问题。但这并不是说明,本文提出的方法仅适用这个最简单的问题。

不稳定性的来源

上述问题的最优解中,

V

a

都是未知的,一般会用估计值

\hat{V}

\hat{a}

。正是这些估计值会导致结果的不稳定性,他们细微的变化会极大的导致结果变化。这种不稳定性可以充以下两个方面说明。

噪音:假设一个TxN的矩阵X,由N的独立同分布的随机变量组成,它们的期望为0,方差为

\sigma^2

。矩阵

C=T^{-1}X^\prime X

有特征值

\lambda

,根据Marcenko–Pastur理论(该定理解释了独立同分布随机变量协方差矩阵的特征值分布情况,这些特征值反映的是各种噪音的波动性),当

N\rightarrow +\infty

T\rightarrow +\infty

1<T/N< +\infty

时,

\lambda

的概率密度函数为:

f[\lambda]=\left\{\begin{array}{ll} \frac{T}{N} \frac{\sqrt{\left(\lambda_{+}-\lambda\right)\left(\lambda-\lambda_{-}\right)}}{2 \pi \lambda \sigma^{2}} & \text { if } \lambda \in\left[\lambda_{-}, \lambda_{+}\right] \\ 0 & \text { if } \lambda \notin\left[\lambda_{-}, \lambda_{+}\right] \end{array}\right.

其中,

\lambda

的最大取值

\lambda_+ = \sigma^2\left( 1+\sqrt{\frac{N}{T}} \right)^2

,最小值

\lambda_- = \sigma^2\left( 1-\sqrt{\frac{N}{T}} \right)^2

。当

\sigma^2=1

时,

C

X

的相关系数矩阵。

但是,实际情况中

\frac{N}{T} \rightarrow 1

,这时

\lambda_-

趋近0,这就导致

\hat{V}

的行列式接近0,

\hat{V}

的逆矩阵就不能很稳健的计算,那么由此得到的解就不稳定。

信号:当相关矩阵为单位矩阵时,特征值函数为一条水平线。 除了这种理想情况下,至少有一个变量子集显示出比其他变量子集更大的相关性,从而在相关矩阵中形成一个簇。当k个变量形成一个集群时,它们更容易暴露于一个共同的特征向量,这意味着相关的特征值解释了更大的数量的方差。但是由于相关性矩阵的迹恰好是N,这意味着一个特征值只能以牺牲该簇中其他K - 1个特征值为代价而增加,从而导致条件数大于1。对于相关性矩阵聚类的特性带来的不稳定性,作者提出了嵌套聚类优化(NCO)

蒙特卡罗模拟法MCOS

MCOS求解w的过程一共包含了五个步骤:

1、估计均值和方差:以 为参数生成矩阵 ,计算矩阵 的均值和方差

参考如下代码:

import numpy as np,pandas as pd
from sklearn.covariance import LedoitWolf

def simCovMu(mu0,cov0,nObs,shrink=False):
      x=np.random.multivariate_normal(mu0.flatten(),cov0,size=nObs)
      mu1=x.mean(axis=0).reshape(-1,1)
      if shrink:cov1=LedoitWolf().fit(x).covariance_
     else:cov1=np.cov(x,rowvar=0)
     return mu1,cov1

2、去噪音:这一步为了解决上文提到的由于噪音带来的不稳定性。首先用KDE算法,将特征值进行Marcenko-Pastur 分布拟合。这样就能够将噪音相关的特征值从信号相关的特征值分离出来。

参考以下代码:

from sklearn.neighbors.kde import KernelDensity
from scipy.optimize import minimize

def fitKDE(obs, bWidth=.25, kernel='gaussian', x=None):
      # Fit kernel to a series of obs, and derive the prob of obs
      # x is the array of values on which the fit KDE will be evaluated
      if len(obs.shape)==1: obs=obs.reshape(-1,1)
      kde=KernelDensity(kernel=kernel, bandwidth=bWidth).fit(obs)
      if x is None: x=np.unique(obs).reshape(-1,1)
      if len(x.shape)==1:  x=x.reshape(-1,1)
      logProb=kde.score_samples(x) # log(density)
      pdf=pd.Series(np.exp(logProb), index=x.flatten())
      return pdf
#------------------------------------------------------------------------------
def mpPDF(var,q,pts):
      # Marcenko-Pastur pdf
      # q=T/N
      eMin,eMax=var*(1-(1./q)**.5)**2,var*(1+(1./q)**.5)**2
      eVal=np.linspace(eMin,eMax,pts)
      pdf=q/(2*np.pi*var*eVal)*((eMax-eVal)*(eVal-eMin))**.5
      pdf=pd.Series(pdf,index=eVal)
      return pdf
#------------------------------------------------------------------------------
def errPDFs(var,eVal,q,bWidth,pts=1000):
      # Fit error
      pdf0=mpPDF(var,q,pts) # theoretical pdf
      pdf1=fitKDE(eVal,bWidth,x=pdf0.index.values) # empirical pdf
      sse=np.sum((pdf1-pdf0)**2)
      return sse
#------------------------------------------------------------------------------
def findMaxEval(eVal,q,bWidth):
      # Find max random eVal by fitting Marcenko's dist to the empirical one
      out=minimize(lambda *x:errPDFs(*x),.5,args=(eVal,q,bWidth),
      bounds=((1E-5,1-1E-5),))
      if out['success']: var=out['x'][0]
      else: var=1
      eMax=var*(1+(1./q)**.5)**2
      return eMax,var
      
#------------------------------------------------------------------------------
def corr2cov(corr,std):
      cov=corr*np.outer(std,std)
      return cov
#------------------------------------------------------------------------------
def cov2corr(cov):
      # Derive the correlation matrix from a covariance matrix
      std=np.sqrt(np.diag(cov))
      corr=cov/np.outer(std,std)
      corr[corr<-1],corr[corr>1]=-1,1 # numerical error
      return corr
#------------------------------------------------------------------------------
def getPCA(matrix):
      # Get eVal,eVec from a Hermitian matrix
      eVal,eVec=np.linalg.eigh(matrix)
      indices=eVal.argsort()[::-1] # arguments for sorting eVal desc
      eVal,eVec=eVal[indices],eVec[:,indices]
      eVal=np.diagflat(eVal)
      return eVal,eVec
#------------------------------------------------------------------------------
def denoisedCorr(eVal,eVec,nFacts):
      # Remove noise from corr by fixing random eigenvalues
      eVal_=np.diag(eVal).copy()
      eVal_[nFacts:]=eVal_[nFacts:].sum()/float(eVal_.shape[0]-nFacts)
      eVal_=np.diag(eVal_)
      corr1=np.dot(eVec,eVal_).dot(eVec.T)
      corr1=cov2corr(corr1)
      return corr1
#------------------------------------------------------------------------------
def deNoiseCov(cov0,q,bWidth):
      corr0=cov2corr(cov0)
      eVal0,eVec0=getPCA(corr0)
      eMax0,var0=findMaxEval(np.diag(eVal0),q,bWidth)
      nFacts0=eVal0.shape[0]-np.diag(eVal0)[::-1].searchsorted(eMax0)
      corr1=denoisedCorr(eVal0,eVec0,nFacts0)
      cov1=corr2cov(corr1,np.diag(cov0)**.5)
      return cov1

3、最优化:根据各种方法计算最优权重,比如CVO或者上文提到的NCO,NCO的代码如下。NCO的方法能够控制信号带来的不稳定性,具体步骤如下:

  • 利用相关性矩阵对变量进行聚类;
  • 对每个子簇进行最优权重计算,这样可以把每个子簇看成一个变量,各子簇之间的协方差矩阵称为简化版协方差矩阵(Reduced Covariance Matrix);
  • 计算各子簇之间的最优权重;
  • 结合上述两个步骤就可以得出每个变量最终的最优权重。

详细代码如下:

from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_samples
#------------------------------------------------------------------------------
def clusterKMeansBase(corr0,maxNumClusters=None,n_init=10):
    dist, silh=((1-corr0.fillna(0))/2.)**.5,pd.Series() # distance matrix
    if maxNumClusters is None:
        maxNumClusters=corr0.shape[0]/2
    for init in range(n_init):
        for i in xrange(2,maxNumClusters+1): # find optimal num clusters
            kmeans_=KMeans(n_clusters=i,n_jobs=1,n_init=1)
            kmeans_=kmeans_.fit(dist)
            silh_=silhouette_samples(dist,kmeans_.labels_)
            stat=(silh_.mean()/silh_.std(),silh.mean()/silh.std())
            if np.isnan(stat[1]) or stat[0]>stat[1]:
                silh,kmeans=silh_,kmeans_
    newIdx=np.argsort(kmeans.labels_)
    corr1=corr0.iloc[newIdx] # reorder rows
    corr1=corr1.iloc[:,newIdx] # reorder columns
    clstrs={i:corr0.columns[np.where(kmeans.labels_==i)[0]].tolist() for \
            i in np.unique(kmeans.labels_)} # cluster members
    silh=pd.Series(silh,index=dist.index)
    return corr1,clstrs,silh
#------------------------------------------------------------------------------
def optPort(cov,mu=None):
    inv=np.linalg.inv(cov)
    ones=np.ones(shape=(inv.shape[0],1))
    if mu is None:mu=ones
    w=np.dot(inv,mu)
    w/=np.dot(ones.T,w)
    return w
#------------------------------------------------------------------------------
def optPort_nco(cov,mu=None,maxNumClusters=None):
    cov=pd.DataFrame(cov)
    if mu is not None:mu=pd.Series(mu[:,0])
    corr1=cov2corr(cov)
    corr1,clstrs,_=clusterKMeansBase(corr1,maxNumClusters,n_init=10)
    wIntra=pd.DataFrame(0,index=cov.index,columns=clstrs.keys())
    for i in clstrs:
        cov_=cov.loc[clstrs[i],clstrs[i]].values
        mu_=(None if mu is None else mu.loc[clstrs[i]].values.reshape(-1,1))
        wIntra.loc[clstrs[i],i]=optPort(cov_,mu_).flatten()
    cov_=wIntra.T.dot(np.dot(cov,wIntra)) # reduce covariance matrix
    mu_=(None if mu is None else wIntra.T.dot(mu))
    wInter=pd.Series(optPort(cov_,mu_).flatten(),index=cov_.index)
    nco=wIntra.mul(wInter,axis=1).sum(axis=1).values.reshape(-1,1)
    return nco

4、蒙特卡罗模拟:结合以上所有步骤,进行多次模拟计算,步骤1每次模拟的 ,都计算出对应的最优解

def monteCarlo(mu0,cov0,nObs,nSims,bWidth,minVarPortf,shrink):
      w1=pd.DataFrame(columns=xrange(cov0.shape[0]),
      index=xrange(nSims),dtype=float)
      w1_d=w1.copy(deep=True)
      for i in range(nSims):
            mu1,cov1=simCovMu(mu0,cov0,nObs,shrink)
      if minVarPortf: mu1=None
      if bWidth>0: cov1=deNoiseCov(cov1,nObs*1./cov1.shape[1],bWidth)
      w1.loc[i]=optPort(cov1,mu1).flatten()
      w1_d.loc[i]=optPort_nco(cov1,mu1,int(cov1.shape[0]/2)).flatten()
      return

5、误差评估:把步骤4计算的 与使用原始均值方差 计算出的最优权重 进行比较,计算误差,误差的定义可以是以下定义之一,或其他任何合理的定义:

a. 均值误差:

\left(w^{*}-\hat{w}^{*}\right)^{\prime} \mu

b. 方差误差:

\left(w^{*}-\hat{w}^{*}\right)^{\prime} V\left(w^{*}-\hat{w}^{*}\right)

c. 夏普误差:

\frac{\left(w^{*}-\hat{w}^{*}\right)^{\prime} \mu}{\left(w^{*}-\hat{w}^{*}\right)^{\prime} V\left(w^{*}-\hat{w}^{*}\right)}

现成的工具包

上文给出的代码多以说明性为目的,在真实研究中应用还有所欠缺,Github上有一个开源的完善的针对本片论文的工具包:

https://github.com/enjine-com/mcos

内部实现了多种最优化算法,包括Markowitz Optimization、 Nested Cluster Optimization、Risk Parity及Hierarchical Risk Parity。请看下面示例说明,针对近20只美股,对不同的权重优化算法进行比较,作者首先使用的ExpectedOutcomeErrorEstimator就是我们上文步骤5提到均值误差评估器。

import numpy as np
import pandas as pd
from mcos import optimizer
from mcos import observation_simulator
from mcos import mcos
from mcos.error_estimator import ExpectedOutcomeErrorEstimator, SharpeRatioErrorEstimator, \
    VarianceErrorEstimator
from mcos.covariance_transformer import DeNoiserCovarianceTransformer, AbstractCovarianceTransformer
from mcos.observation_simulator import AbstractObservationSimulator, MuCovLedoitWolfObservationSimulator, \
MuCovObservationSimulator
from pypfopt.expected_returns import mean_historical_return
from pypfopt.risk_models import sample_cov
import warnings
warnings.filterwarnings('ignore')
# Create dataframe of price history to use for expected returns and covariance
def prices_df() -> pd.DataFrame:
    tickers = ['goog','baba', 'amzn', 'wmt', 'glpi', 'bac', 'uaa', 'shld', 'jpm', 'sbux', 'amd', 'aapl','bby',
              'ge', 'rrc', 'ma','fb']
    total_df = pd.DataFrame()
    for id in tickers:
        temp = pd.read_csv( id + '.us.txt', parse_dates=True, index_col='Date')
        temp = pd.DataFrame(temp['Close']).rename(columns={"Close":id})
        if total_df.empty:
            total_df = temp
        else:
            total_df = total_df.join(temp)
    return total_df
# Choose the number of simulations to run
num_sims = 50
# Select the optimizers that you would like to compare
op = [optimizer.HRPOptimizer(), optimizer.MarkowitzOptimizer(),optimizer.NCOOptimizer(), optimizer.RiskParityOptimizer()]
# select the metric to use for comparison
ee = ExpectedOutcomeErrorEstimator()
# select your optional covariance transformer
cov_trans = DeNoiserCovarianceTransformer()
# convert price history to expected returns and covariance matrix
mu = mean_historical_return(prices_df()).values
cov = sample_cov(prices_df()).values
# select your observational simulator
obs_sim = MuCovObservationSimulator(mu, cov, num_sims)
# Run the simulation
results = mcos.simulate_optimizations(obs_sim, num_sims, op, ee, [cov_trans])
print(results)
results.plot.bar()

上图为利用均值误差评估器,对各权重优化模型评估的结果,我们可以发现Risk Parity模型表现得最稳健。

该工具包还支持自定义优化器,并对其进行评估,有兴趣的Quant可以尝试~

量化投资与机器学习微信公众号,是业内垂直于量化投资、对冲基金、Fintech、人工智能、大数据等领域的主流自媒体。公众号拥有来自公募、私募、券商、期货、银行、保险、高校等行业20W+关注者,连续2年被腾讯云+社区评选为“年度最佳作者”。

本文分享自微信公众号 - 量化投资与机器学习(Lhtz_Jqxx),作者:全网Quant都在看

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2021-07-15

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 基石长固,腾云共舞——腾讯第二届云+数据中心分享日纪实

    去年九月,腾讯数据中心在北京腾讯汇举办了第一届数据中心分享日。与会的有互联网公司BAT的代表、工信部电信研究院、三大运营商以及腾讯的战略合作供应商等。在第一届分...

    腾讯数据中心
  • TensorFlow v2.x使用说明[1]-概要与更新

    版权声明:本文为博主原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明。 ...

    锦小年
  • 太厉害了!Seaborn也能做多种回归分析,统统只需一行代码

    lmplot是一种集合基础绘图与基于数据建立回归模型的绘图方法。通过lmplot我们可以直观地总览数据的内在关系。显示每个数据集的线性回归结果,xy变量,利用'...

    数据STUDIO
  • 对于小批量梯度下降以及如何配置批量大小的入门级介绍

    随机梯度下降是训练深度学习模型的主要方法。

    用户1284604
  • 互联网医疗健康产业联盟酝酿大动作 这八份礼物一定有你要的

    从2014年“春江水暖”,到2015年成为投资界“香饽饽”,到2016年“健康中国”成为国家战略,如今,互联网医疗健康产业正在逐步回归实质——借助互联网手段,实...

    企鹅号小编
  • 通过深度学习识别和验证基于脑额叶区-后叶区功能失衡的重大精神疾病内的亚型

    精神分裂症(SZ)、双相情感障碍(BD)和重性抑郁症(MDD)是在精神疾病领域常见的三种疾病,合称为重大精神疾病(MPD),长期以来都是依据不同的核心症状被作为...

    悦影科技
  • 通过深度学习识别和验证基于脑额叶区-后叶区功能失衡的重大精神疾病内的亚型

    精神分裂症(SZ)、双相情感障碍(BD)和重性抑郁症(MDD)是在精神疾病领域常见的三种疾病,合称为重大精神疾病(MPD),长期以来都是依据不同的核心症状被作为...

    悦影科技
  • 高三学生发表AI论文,提出针对网络暴力问题的新模型AdaGCN

    【导读】近日,在清华大学举行的丘成桐中学科学奖半决赛落下帷幕,来自海内外的 72 支队伍获得了总决赛的入场券,北京师范大学附属实验中学的高三学生白行健,也在其中...

    AI科技大本营
  • 一步步靠近:Rust入门小百科

    从2015年Rust 发布1.0 版本以来,Rust 语言已经被广泛应用于各大公司及诸多领域。每一年,Rust 社区都会聚集在一起制订路线图,规划Rust 未来...

    用户1682855
  • R2015b 版本

    R2015b 版本 MATLAB 产品系列更新: MATLAB: 新增更快运行 MATLAB® 代码的执行引擎;用于创建、分析图形和网络并实现可视化的图形...

    万木逢春
  • 美年健康复牌点评:美年健康业绩高速增长,携手慈铭实现强强联合

    携手慈铭,健康体检领域两大龙头公司强强联合:公司通过三步来完成对慈铭体检的收购,交易价格公允;携手慈铭,实现网点互补,运营层面以期实现规模效应,经营层面以期解决...

    曾响铃
  • 可以丢掉SGD和Adam了,新的深度学习优化器Ranger:RAdam + LookAhead强强结合

    给大家介绍一个新的深度学习优化器,Ranger,同时具备RAdam和LookAhead的优点,一行代码提升你的模型能力。

    AI算法与图像处理
  • 算法--基础

    版权声明:本文为博主原创文章,遵循 CC 4.0 by-sa 版权协议,转载请附上原文出处链接和本声明。

    奋飛
  • 架构思想

    Java高级架构
  • 架构思想

    关于什么是架构,一种比较通俗的说法是 “最高层次的规划,难以改变的决定”,这些规划和决定奠定了事物未来发展的方向和最终的蓝图。 从这个意义上说,人生规划也是一种...

    Java高级架构
  • 人工智能普及,你要失业了?

    IT大咖说
  • 资产瞎配模型(二):对瞎配(一)中净值计算错误的纠正

    上上周发的那篇资产瞎配模型,事实证明,果然是瞎配,有大佬指出组合净值计算有一定的问题,所以这里对净值计算部分及进行改正,重新计算结果。

    量化小白
  • 波动率目标策略,没有想象的那么简单!

    量化投资与机器学习公众号独家解读 量化投资与机器学公众号 QIML Insight——深度研读系列 是公众号今年全力打造的一档深度、前沿、高水准栏目。

    量化投资与机器学习微信公众号
  • 2019医疗大数据企业排行榜

    近日,互联网周刊发布“2019医疗大数据企业排行榜”,一起来看一下哪些公司上榜了?

    华章科技

扫码关注云+社区

领取腾讯云代金券