风险模型在时间序列上的改进 ——《系列之三十一》

      假设是一个p维随机向量，有n个观察值，随机向量的协方差记为。如果满足p维正态分布，则 Σ 的极大似然估计正好等于样本协方差矩阵：

S 是一个无偏估计，如果 p 是一个不变的常数，则在时， 。不过在估算股票收益率协方差矩阵时，由于股票收益率的分布在时间序列方向上不稳定，不建议用太长的历史数据，所以经常碰到的情况是用过去一年252个交易日的数据去估算1000甚至更多只股票的协方差矩阵，样本协方差矩阵S的秩等于 n – 1，因此当p≥n时，S不可逆；即使 p <n,根据Ledoit(2004)里面的引理2.1，S特征值也会比 Σ 特征值的分布区间广，因此S特征值的最大值比Σ 特征值的最大值大，S特征值的最小值比Σ 特征值的最小值小，造成矩阵S的条件数（Condition Number, 矩阵最大特征值除以最小特征值）偏高，偏高幅度与p/n 的大小正相关。后续的组合优化过程中可能会对S求逆或者以S为系数矩阵求解线性方程，S条件数过大会让计算结果对数据误差十分敏感，使得组合优化得到的组合权重不稳定。

        解决这个问题的常用方法有三种：
        1） 增加样本数量n。前文提到不宜用太久之前的历史数据，因此可以考虑提高数据的采样频率，也就是采用高频数据。不过高频数据在增加样本数量的同时也在增加数据噪音，统计工具上可能要有所改变。另外高频数据计算得到的是高频协方差矩阵，高频数据在收益率序列自相关性和低频数据可能有差异，如何把高频协方差矩阵转换成低频协方差矩阵也是需要考虑的问题。
        2） 降低变量维度p。BARRA因子模型属于这一类，也称作结构化风险模型。它通过寻找少量因子来解释股票收益率，把估计股票收益率协方差矩阵转换成估计因子收益率协方差矩阵和个股特质方差，把需要估计的参数数量从大幅降为O(p)，减小估计量的方差。因子对股票收益率的解释度越高，降维效果越明显。
        3） 通过函数转换调整样本协方差矩阵的特征值，同时提高估计准确性。这也是下文要介绍的压缩估计方法。

2.2 线性与非线性压缩估计

        首先介绍压缩估计量方法，这个在后面的因子模型里也会用到。对于p维对称矩阵和，常用Frobenius 距离来衡量两个矩阵的接近程度：

      如果Σ 是要估计的真实协方差矩阵，是某个协方差矩阵估计量，可以定义估计量的估计误差为

        样本协方差矩阵S 的条件数偏大，也就是说S特征值的最大值和最小值距离太远，需要压缩来降低条件数。一种简便的线性压缩方法是给S加上一个目标矩阵T，目标矩阵的选取有很多种（Schafer 2005），以单位阵I为例，线性压缩（LS，Linear Shrinkage）估计量可以表示成：

线性压缩估计即是要在所有可能的线性组合中，找到一个估计误差最小的，即求解

可显式求得这个问题的解为

其中μ等于Σ 特征值的平均值，表示样本协方差矩阵的估计误差，度量的是样本协方差矩阵特征值分布的离散程度。Σ 的真实值不知道，也就无法计算，不过在一定的理论假设下，Ledoit(2004)构建了的样本统计量，保证在时，对应的样本压缩估计量收敛于。
        样本协方差矩阵是（1）式的一个特例（ρ=0），因此LS估计量的估计误差小于等于样本协方差（n,p足够大的时候），误差减小的幅度取决于两个数值p/n 和 ρ，值越大，估计误差减小幅度越明显。LS估计量不依赖于随机变量的分布，而且中间过程涉及的变量都可以显式计算，因此计算效率非常之高，全市场3000只股票计算一次协方差矩阵耗时在毫秒级别，因此也是我们之前研究中最常用的协方差矩阵估计方法。Python（sklearn. covariance.LedoitWolf）, MATLAB（作者个人主页有作者自己编写的代码可下载），R（nlshrink package）都有成熟的工具包可调用。
        记S的特征值为，其对应特征向量为，则S可以谱分解为。压缩估计量可以表示成：

可以看到，S的特征值被线性压缩到了的，也就是说向真实协方差矩阵Σ 特征值的平均值靠拢，靠拢程度取决于ρ的大小。这样特征值的分布区间相对S而言得到了一定的压缩，矩阵的条件数减小，用它做后续的组合优化，结果对数据误差的敏感性也会降低。
        LS是给样本协方差矩阵S所有特征值设定了同一个目标，做全局压缩。另一种更灵活的非线性压缩（NLS, Nonlinear Shrinkage）方式由Ledoit(2017)近期提出，它的形式和(2)很类似，可以写作:

这里d(⋅) 是一个一元函数，在LS的（2）式中d(⋅)是一个一元函数。在对总体的协方差矩阵，数据产生过程和函数d(⋅)的特性做了一些理论假设后，可以证明下面oracle 函数可以最小化NLS估计量的估计误差：

其中 c=p/n, f(⋅)是协方差矩阵特征值的极限谱分布函数, 是其对应的Hilbert变换。Ledoit(2017)通过核函方法计算得到样本协方差矩阵谱分布函数、对应的Hilbert变换和oracle函数，并证明当时。
        NLS和LS对样本协方差矩阵特征值的压缩效果不同。对于极大和极小的特征值，NLS和LS一样也会将其往均值方向压缩，但中间区域的特征值压缩方向有可能是远离均值的，NLS会把特征值向其附近特征值取值集中的区域压缩，每个不同特征值处的压缩强度都是不一样的，是一种局部性质。Ledoit (2017)通过仿真模拟发现，在真实协方差矩阵的条件数较大、样本集中度（p/n）较小或样本数量较多时，NLS相对LS的估计误差减小幅度非常明显。而且NLS计算过程里涉及的变量也都是可以显式计算的，运算速度非常之快，和LS基本相当，因此之前采用LS估计的研究人员可以非常便捷的转移到NLS上。
        需要注意的是，Ledoit(2017) 附录 C有一处小错误，其文后所附的MATLAB代码需要做几处对应的修改。

2.3 因子模型
        假设市场上有p个股票，我们已经找到K个风险因子，股票 i 在因子k上的暴露度为, 记p×K矩阵B为因子暴露度矩阵，则做横截面回归

可以估算得到横截面上的风险因子纯因子收益率。考虑到横截面上股票特质方差的差异，需要采用WLS（Weighted Least Square）以保证得到的估计量同时也是极大似然估计，权重为个股特质方差的倒数。个股特质方差和其市值平方根倒数近似成正比（参考BARRA USE4和我们之前报告的A股实证结果），所以实际计算的时候可以采用个股市值平方根作为权重。由此，股票的协方差矩阵可如下计算：

每个横截面上做一次回归得到f和ϵ的时间序列后，F和S可以在时间序列方向上进行估计。
        横截面回归的Rsquared 可以用来度量风险因子对股票收益的解释程度。不过要注意的是，Rsquared的计算方式和测试股票池不一样，得到的数值会有很大差别。由于横截面回归使用了WLS方法，因此可以计算weighted rsquared，不过这个数值会比传统OLS里面的adjusted rsqaured高很多，例如我们参考CNE5 构建的因子模型，全市场回归得到的月度weighted rsqaured 平均超过40%，但 adjusted rsquared 只有20%左右。因子模型在大股票里的解释能力更强，如果缩减股票池，weighted rsquared 可以超过50%，adjusted rsquared 可以超过30%。在协方差矩阵估计中，adjusted rsquared 的参考意义更大，因为协方差矩阵估计过程中，各个股票的数据地位是平等的，adjusted rsquared 越大，能被风险因子解释的方差越多，因子模型的降维效果越明显。
        我们在去年的专题报告《A股市场分析中》参考BARRA CNE5文档构建了一个因子风险模型，并和LS做对比，当时实证发现因子风险模型的效果不如LS。不过严格的讲，当时报告里的实证方法不公平，因子模型采用的是月频的风险因子数据，LS模型用的则是日频收益率数据，为保证公平性，报告后文的实证中，风险因子数据都是采用日频更新，因子收益率也是每日计算。
        在上述标准模型基础上，我们用非线性压缩方法来估计因子收益率协方差矩阵F和特质方差矩阵S，以降低样本协方差矩阵条件数过大的问题， CNE5中对特质方差矩阵的Bayesian Shrinkage是一种线性压缩估计方法，采用非线性压缩有可能进一步提升模型表现。另外我们给F和S也加上了一个EWMA架构，通过调整衰减系数来调节模型对近期数据的敏感性。

2.4 多元GARCH类模型
        压缩估计和标准的因子模型相对样本协方差矩阵的改进都属于横截面方向上的，也就是说它们都基于股票收益率在时间序列方向上是独立同分布的假设，不过股票在时间序列上的异质性明显，可以考虑从这个方向上改进风险模型。在一维情形下，我们可以用GARCH类模型来描述条件方差的动态变化，因此一个最直接的思路是将GARCH模型推向高维。
        记t 时刻 p 个股票的收益率为 p 维向量 ，它基于t-1时刻已有市场信息的条件预期收益和协方差分别假设为。我们需要赋予一定的结构来计算极大似然函数，目前最常用的结构有两种：Engle(1995)的BEKK模型和Engle(2002)的DCC（Dynamic Conditional Correlation）模型。
        BEKK模型直接对协方差矩阵建模，假设满足如下递归形式：

        DCC则是对相关系数矩阵建模，假设是一个p×p 对角阵, 其第i个对角线元素是股票i 在t时刻的条件方差，DCC模型假设相关系数矩阵满足如下形式：

其中。
        在满足正态分布的假设下，容易计算得到对数极大似然函数为

对于低维度问题，p比较小，上述极大似然函数比较好计算，R里面有成熟的package（rmgarch）可以直接调用。但如果p较大，每次求的逆需要耗费运算，极大似然函数的求值将十分费时，而且后面用数值方法优化极大似然函数时，极大似然函数可能会被求值成百上千次，个人电脑上可能需要几天的时间才能算完，实用性有限。
        Pakel(2017) 借鉴Lindsay(1988)提出的CL（Composite Likelihoods）思想大大降低了DCC和BEKK模型估计参数的计算量，并在一定理论假设前提下，保证估计量的一致性。这个方法分为两部，第一步是用传统样本协方差方法来估计BEKK模型里的Σ 和DCC里面的C，第二部计算CL函数，和(3)式的极大似然函数相比，它把p个随机变量两两配对，计算二维的对数极大似然函数再进行累加。这样矩阵的逆运算都是二维的，p较大时，运算速度会快很多倍。Engle(2017)对第一步的方法做了改进，采用NLS方法估计Σ和C，这样做的好处是在股票数量较多时，可以大幅降低估计误差。
        我们用Engle(2017)的方法编写python代码在A股进行了测试，时间点是2009.01.04，用过去一年的日频收益率数据分别在沪深300成份股（HS300）、中证500成份股（ZZ500）和小市值500（SC500, 中证800之后市值最大的500只股票）进行了参数估计，结果如图1所示。BEKK模型的alpha参数大概0.02附近，beta参数在0.98附近，这个和Pakel(2017)的在标普500里面的参数估计结果很接近，他们估算得到alpha大概在0.03附近，beta 在0.96附近。DCC模型估计得到的alpha和美股相差较大：A股估计得到的数值大概在0.001左右，而美股大概在0.007，也就是说A股股票间的条件相关系数矩阵对新收益率数据不敏感，动态性不强；造成这种结果的一种可能性是A股股票间的相关系数确实变化不大（这个与主观感觉有出入），另一种可能性是DCC模型设定不适合A股，存在模型设定偏误。 DCC模型的beta在美股和A股都在0.99左右，比较接近。

模型参数估计耗费的时间值得关注。在个人电脑上（i5-4590+16G内存+pycharm+  python3.5），估算沪深300成份股300只股票的DCC和BEKK模型参数，只需10分钟左右，但增加到500只股票（ZZ500&SC500）时，耗时接近一个小时。CL方法和最初的MLE估计方法几天的耗时相比已经大幅降低了，而且在实际实盘投资中也可用，开盘前运行一次，输入到组合优化中得到组合权重即可。不过对于长时间的历史回溯而言，累加起来的耗时量非常惊人，调参数困难。程序耗时的原因可能和代码编写有关系，计算CL函数时需要通过BEKK和DCC的迭代式反复循环计算，而后续在优化CL函数时，CL函数可能会被反复求值上千次，总运算量巨大。虽然迭代式的代码我们已经用Cython进行了加速，但目前看来效果并不理想，用纯C语言来写可能效果会更佳。
       鉴于DCC和BEKK模型的耗时，我们可以考虑采用DCC模型的一个简化版本CCC GARCH，它假设DCC模型中的条件相关系数矩阵是一个常数矩阵，不随时间变化，条件协方差的变化完全由个股条件方差的变化引起。这样估计参数时，只需要对每个股票做时间序列上的一维GARCH模型拟合，然后计算devolatised return的相关系数矩阵即可，计算量大幅降低。图1在A股估计得到的DCC模型alpha很小，在时间序列上的动态性不强，支持我们采用CCC GARCH模型做简化。Pakel(2017)做过大量Monte-Carlo仿真测试，结果显示当股票数量较少时（100只左右），DCC相对CCC模型的估计误差降幅明显，但随着股票数量的逐步增加（1000只），误差的相对降幅越来越小；而且如果用真实的股票收益率数据构建全局最小方差组合（GMVP），两者的结果基本相当。
        最后补充一点，Pakel(2017)和Engle(2017)在做实证时都是选取历史上一直都在正常交易的股票作为标的，但A股实际投资中，不可避免会碰到很多历史长期停牌的股票，我们研究发现缺失数据的填补方式对模型的使用效果影响较大。