基于 DolphinDB 的多因子风险模型实践

/* winsorized
对因子表进行winsorized处理
因子表必须包含三列：记录日期、股票代码和原始因子，其中原始因子必须是第三列
Input :
         tbName
Output:  
         factor table after winsorized */

MAD 方法的具体实现步骤如下：

1. 计算数据集的中位数（Median）作为数据的中心位置。

2. 对每个数据点，计算其与中位数的绝对偏差（Absolute Deviation）。

3. 计算所有绝对偏差的中位数，即 Median Absolute Deviation。
   
   
4. 确定离群值的阈值。通常，离群值被定义为与中位数的偏差超过一定倍数（如 3 倍）的 MAD 的数据点。
   
   
5. 鉴定离群值。对于超过阈值的数据点，可以将其标记为离群值或进行进一步的分析。
   
   

市值加权标准化是为了去除极端市值对因子计算贡献产生影响（CNLT ）。该方法是为了将因子量纲统一在相接近的水平。对于第 n 个个股的第 k 个原始因子暴露 x_n,k^Raw，标准化流程为：

其中 μ_k 和 σ_k 分别是对应因子的市值加权均值、等权标准差。(CNLT)

本文市值加权标准化对应接口函数 standardized，其中 adjusted=true 表示采纳市值中性化处理。

/* standardized
对因子表进行标准化和市值中性化处理
因子表必须有三列：记录日期、股票代码和原始因子，原始因子必须是第三列
Input :
         tbName
         adjusted    Market-neutralize or not
Output:  
         factor table after standardized */

以上两种接口函数 winsorized, standardized 嵌入在 getAllFactors 函数中，getAllFactors 通过指定参数 normlizing 和 scaling 分别实现标准化和去极值的过程。此外在得到行业因子、风格因子后，基于 getAllFactors 函数可以得到合并了上述所有因子的宽表，可进一步用于单因子的有效性检验。

/* getAllFactors
获取所有因子
Input:   normlizing           true   (Default)   标准化
         scaling              true   (Default)   去极值
         decap                true   (Default)   市值中性化
         industry_weighted    true   (Default)   行业市值权重加权
         industry_method     'CITIC' (Default)、'SW_2021'
         startTime         2022.01.03(Default)
         endTime           2023.01.02(Default)
Output:  
         factor table */

Factors = getAllFactors(st=st,et =et, normlizing = true,
                        scaling = true,decap = false,industry_method = 'CITIC',
                        industry_weighted = false)
select * from Factors limit 100

/* getRegTable
 获取用于回归的经过处理的回归因子表、包含个股收益率、因子暴露、行业因子、行业变量、回归权重
Input:   
        factorsTable         false  (Default)          是否使用提供的初始因子表
        tbName               NULL   (Default)          初始因子表
        normlizing           true   (Default)          标准化
        scaling              true   (Default)          去极值
        decap                true   (Default)          回归市值中性化
        industry_weighted    true   (Default)          行业因子加权重
        industry_method     'CITIC' (Default)、'SW_2021'
        st                2022.01.03(Default)
        et                2023.01.02(Default)
Output:  
        regression table  
 */

通过调用：

fTable = getRegTable(factorsTable = true,tbName = Factors,st= st,
                     et = et,normlizing = normlizing ,scaling = scaling,
                     decap = decap, industry_method = industry_method, 
                     industry_weighted = industry_weighted)

执行上述命令可以实现如下效果：

处理之前，getAllFactors 得到的原始三级因子宽表：

处理之后，getRegTable 处理后的三级因子宽表：

• 针对风格因子 s，单因子的检验回归模型如下：
  
  

  此处使用（Weighted Least Squares，加权最小二乘）以 r_n,t 为因变量，以风格因子 s 的因子载荷 x_n,s,t 为自变量，预测 f_s.t（风格因子的收益率）。

• 针对行业因子 i，单因子的检验回归模型如下：
  
  

  此处使用以 r_n,t 为因变量，以行业因子 i 的因子载荷 x_n,i,t为自变量，预测 f_i.t（行业因子的收益率）。

  本文的回归模型的接口为 getOneFactorValidate，函数嵌套在 styleValidate 风格因子检验函数、industryValidate 行业因子检验函数中。

  /* getOneFactorValidate
  获取moving wls单因子回归结果统计量的聚合函数
  Input:   y          因变量
           x          自变量
           w          权重
  Output:  
           wls stat  "beta","tstat","R2","AdjustedR2", "Residual" */

T 值检验是检验对应变量的因子收益是否显著区别于 0，可以用来衡量因子的有效性以及一致性。其计算方法为：其中 SE( f ^{^}_i,t) 为估计因子收益 f ^{^}_i,t 的标准误差。本文的 T 检验主要基于 getOneFactorValidate 接口的回归模型得到的 t 统计量进行评估。

在构建多因子模型时，必须要考虑到因子暴露的稳定性。如果因子暴露矩阵每次计算时的变化特别大，那么该模型的稳健性较差。故本文引入因子稳定性系数（Factor Stability Coefficient, FSC）指标，其定义为这个月的因子暴露矩阵与下个月因子暴露矩阵的相关性。

通常情况下，FSC 指标是通过比较从不同数据集或不同时间点上进行的、两个独立的因子分析所得到的因子载荷来进行计算的。本文通过 getFactorsValidation 函数基于斯皮尔曼相关系数（Spearman's Rank Correlation Coefficient）计算 FSC 指标，令 d_t 为因子暴露 x_n,i,t 的第 t 期与第 t+1 期的位次值之差。，则 FSC系数计算方法为：

研究人员通常使用阈值来解释 Factor Stability Coefficient。例如，FSC 大于 0.8通常被认为表示较好的稳定性或一致性，而低于 0.5 的值可能表示稳定性较差。

IC 值被认为在 0.03 上时与市场存在一致或相反的波动规律，此时该因子存在一定的有效性。

本文的单因子有效性检验的对应接口函数 getFactorsValidation。该函数通过输入计算得到的全因子表，给出了每一个单因子模型的单因子收益，t 统计量，拟合优度，FSC 指标和 IC 值。通过这些值可以全面系统地对所有因子进行分析，如下列图所示：

Input:   
         factorsTable         false  (Default)   是否使用提供的因子表
         tbName               NULL   (Default)   回归因子表
         normlizing           true   (Default)   标准化
         scaling              true   (Default)   去极值
         decap                true   (Default)   市值中性化
         industry_weighted    true   (Default)   行业因子加权重
         industry_method     'CITIC' (Default)、'SW_2021'
         st                   2022.01.03(Default)
         et                   2023.01.02(Default)
Output:  
         factor test table  factor_return、tstat、R2、fsc、IC
         */

• 步骤一：首先获取所有因子的检验指标。

   // 获取单风格因子有效性、一致性、稳定性检验
  factorsValid = getFactorsValidation(factorsTable = true,tbName = out,st=2022.01.03,
                                      et =2023.01.02, normlizing = true,scaling = true,
                                      decap = true,industry_method = 'CITIC',
                                      industry_weighted = true)

• 步骤二：绘制因子的 FSC 月频时序图，评价因子稳定性。

  tmp = select record_date,valueType.regexReplace("_stat","") as valueType,
               fsc from factorsValid 
  tmppivot = select fsc from tmp pivot by record_date,valueType
  tbfsc = sql(select = sqlCol(tmppivot.columnNames()[11:20]),from = tmppivot).eval()
  plot(tbfsc,tmppivot.record_date,extras={multiYAxes: false},title = "因子fsc 月频时序图")

  
  
  

  FSC 大部分处于 0.8 以上的因子被公认为因子具备较高的稳定性。上图中除 em, dastd, dtoa 外的其他因子具备较好的稳定性。

• 步骤三：绘制因子的 IC 月频时序图，评价因子一致性。IC 值被认为在 0.03 上时与市场存在一致或相反的波动规律，此时该因子存在一定的有效性。

  tmp1 = select record_date,valueType.regexReplace("_stat","") as valueType,
                abs(ic) as ic from factorsValid 
  tmppivot1 = select ic from tmp1 pivot by record_date,valueType
  tbic = sql(select = sqlCol(tmppivot1.columnNames()[2:10]),from = tmppivot1).eval()
  baseline = take(0.03,(shape tbic)[0])
  plot(table(tbic,baseline),tmppivot1.record_date,
      extras={multiYAxes: false},title = "因子ic 月频时序图")

  
  
  

  可以观察到 atvr（Annualized Traded Value Ratio，年交易比值）因子在 2017 年前保持着较强的相关性，而 2018-2022 年期间 8 个因子呈现周期性起伏的规律。

• 步骤四：绘制因子 t_stat，评价因子有效性。

  tmp2 = select record_date,valueType.regexReplace("_stat","") as valueType,
                tstat from factorsValid 
  tmppivot2 = select tstat from tmp2 pivot by record_date,valueType
  tbstat = sql(select = sqlCol(tmppivot2.columnNames()[11:20]),from = tmppivot2).eval()
  
  baseline_neg = take(-0.03,(shape tbstat)[0])
  baseline_pos = take(0.03,(shape tbstat)[0])
  plot(table(tbstat,baseline_neg,baseline_pos),tmppivot2.record_date,
      extras={multiYAxes: false},title = "因子t_stat 月频时序图")

  
  
  

  图中大部分因子的 t 值距离基线 0.03 和 -0.03 有相当的距离。例如 dastd 和 cmra，一个显著高于 0.03，另一个则显著低于 -0.03。由此可得结论：在 18 年之前，这两个因子与市场存在着显著的正负相关性。

• 对已有因子数据做预处理，可以调用因子预处理的 winsorized、standardized 等函数以实现对风格因子实现极值化、标准化处理。
• 对上述步骤处理后的因子数据处理缺失值，避免因子缺失值对于回归模型结果的影响。

假设最后的因子表为 Factors 表，则可以通过如下方式调用 getRegTable 函数进行缺失值处理。

fTable= getRegTable(factorsTable = true,tbName = Factors,st= st,et = et)

假设最后的因子表为 Factors 表，则可以通过如下方式调用 getFactorsValidation函数进行单因子模型检验，进一步筛选有效因子。

factorsValid= getFactorsValidation(factorsTable = true,tbName = Factors,st= st,et = et)

在指定与自定义因子一一对应的一级因子、二级因子后，通过getFSLevelFactor函数生成基于自定义因子合成后的一级、二级因子。

getFSLevelFactor(fTable,factorsValid,firstFactors,secondFactors,false , "ir",level = "F")
getFSLevelFactor(fTable,factorsValid,firstFactors,secondFactors,false , "ir",level = "S")

贝叶斯收缩（Bayesian Shrinkage）是一个常见的将先验和样本估计值结合起来的手段。进行贝叶斯收缩调整的主要原因是，基于多因子风险模型使用样本内数据计算出的特异性波动率在样本外的持续性很差，极高或极低的个股波动率可能会出现均值回归的情况，那么波动率会被高估或低估。

本文的贝叶斯收缩的接口为 BayesShrinkage 函数，在 getRetTable 函数中当参数 shrink=true 时被调用。

/* BayesShrinkage
BayesShrinkage调整特质波动率
Input:   
         cov                    特质收益率协方差阵
         weight                 市值权重
         q                      λ_F 压缩系数
         
Output:  
         cov                    BayesShrinkage调整后的协方差阵
*/

上述接口的实现步骤如下：

• 步骤一，计算样本估计值——基于 WLS 预测得到的特异性收益的风险 σ^{^}_n。
• 步骤二，计算先验值——首先选择与目标个股处在同一市值的所有股票，再基于市值加权计算的特异性收益率的波动的平均值。
  
  

  其中多因子风险模型把所有个股按照市值分成十档，s_n 表示对应的市值档位，σ^{^}_n 为同处在 s_n 市值档位下的 n 支股票的特异性收益率的波动，w_n 表示 n 支股票按照市值计算的权重。

• 步骤三，贝叶斯收缩（使样本估计值向先验靠拢），计算后验估计。
  
  

  v_n 是在收缩时赋予先验的权重（称为收缩强度系数），具体公式如下：

  
  
  

  其中 N(s_n) 是当前个股市值档位中的个股总数，为样本特异性波动率与先验的偏离程度，为 s_n 市值档位下样本特异性波动率与先验的偏离程度的标准差，q 为经验压缩系数，其意义在于更偏重组内的标准差还是更偏重股票相对组内均值的偏差。当越高，即目标股票的特异性收益率的波动的偏离程度越高，此时样本估计值越不靠谱，因此赋予先验的权重 v_n 则越高。

假设投资组合中的每只股票的特质因子收益率与共同因子收益率不相关，并且每只股票的特质因子收益率也不相关，因此若 w^T表示投资组合中各股票的权重，则投资组合 P 的风险矩阵可表示为如下表达式：

• 行业中性

  行业中性是指多头组合的行业配置与对冲基准的行业配置相一致。行业中性配置的目的在于剔除行业因子对策略收益的影响，仅考察行业内部个股的超额收益。假设 H 为样本股票的行业因子哑变量矩阵，h 为沪深 300 中 30 个行业的对应权重，那么行业中性的权重 w 满足：

  
  
  
• 风格因子中性

  风格因子中性是指多头组合的风格因子较之对冲基准的风险暴露为 0。风格因子中性配置的目的是消除投资组合与市场风格因子之间的风险暴露，使投资组合的收益主要来源于阿尔法收益，而不是某一特定的市场风格。假设 X 为样本第 k 个因子的载荷截面， w_bench 为沪深 300 指数对应权重，那么因子 k 的风格中性权重 w 满足：

  
  
  
• 现金中性

  现金中性是指多空市值保持一致，不留方向性敞口。则现金中性的权重 w 满足：

  
  
  
• 最小收益预测

  控制最小预测收益是多因子模型中组合权重优化的一个约束条件，其目的是确保投资组合具备一定的最低收益水平，并帮助控制风险、避免非理性配置以及维持模型一致性。这样可以使投资组合更符合投资目标和风险偏好，并提高整体投资组合的稳定性和可预测性。假设 r_min为基准收益率，则投资组合的权重 w 应满足：

  
  
  

本文给出行业基准中性和风格基准中性下的控制最小预测收益并最小组合风险方法，对应接口函数 getOptimizeWeights，即实现如下的优化目标。在此接口函数中，可以通过 deIndustry 和 deStyle 指定是否行业中性和风格中性：

/* getOptimizeWeights
组合权重优化中的聚合函数
Input:   
         covf              因子收益协方差矩阵
         delta             特质收益协方差矩阵
         st                2022.01.03(Default)                
         et                2023.01.02(Default)
         ret               预期收益
         r                 0.05                               设置的最小预测收益
         tbName            因子暴露
         deIndustry        true                               行业中性
         deStyle           true                               风格中性
 
Output:  
         weightTable       返回预测的资产配置权重
*/

基于 getOptimizeWeights 接口，可以通过如下脚本实现在控制最小收益下的最小化风险：

optionCode = exec stock_code from getPredicOut(facTable1).predict_ret 
             order by return_day desc limit 20           
// 初步筛选stock1
optionCode = exec stock_code from getPredicOut(facTable2).predict_ret 
             order by return_day desc limit 20  

// 控制收益、最小化风险模型
portWeight1 = getOptimizeWeights(facTable = facTable1,retOut = retOut1,
                                 st = st,et = et, method ="minRiskControlRet",
                                 r = 0.05,optionCode = optionCode)      
// 获得权重组合
portWeight2 = getOptimizeWeights(facTable = facTable2,retOut = retOut2,
                                 st = st,et = et, method ="minRiskControlRet",
                                 r = 0.05, optionCode = optionCode)

index_code = '000300'
CodePre = set(exec stock_code from getPredicOut(facTable1).predict_ret 
              order by return_day desc limit 200)          
 // 初步筛选stock2
CodeWeight = set(exec stock_code 
                from getBenchMark(st=st,et=et,code = index_code) 
                where i_weight != 0)
CodeFac =set(exec stock_code from facTable1 )
optionCode = (CodePre&CodeWeight&CodeFac).keys()
portWeight3 = getOptimizeWeights(facTable = facTable1,retOut = retOut1,
                                 st = st,et = et, method ="minRiskControlRet",
                                 r = 0.005,deStyle = true,optionCode = optionCode)  
// 获得权重组合,并实现在风格上的风险敞口为0
portWeight3 = getOptimizeWeights(facTable = facTable1,retOut = retOut1,st = st,
                                 et = et, method ="minRiskControlRet",r = 0.005,
                                 deIndustry = true,optionCode = optionCode)  
// 获得权重组合,并实现在行业上的风险敞口为0
portWeight4 = getOptimizeWeights(facTable = facTable2,retOut = retOut2,st = st,
                                 et = et, method ="minRiskControlRet",r = 0.05,
                                 optionCode = optionCode)

事后资产配置指在实际收益数据可用之后，根据实际的历史收益数据进行的资产配置。这个过程发生在投资决策之后，基于实际观察到的历史收益数据对资产进行重新配置。因此根据市值或者是等权法评估已有指数的 Bias，可以计算出指定组合的偏差统计量和 Q 统计量，以对事后资产配置进行评估。

• 因子（组合）Bias
• 资产（组合）Bias

本文基于 getFacSpecialBias 函数，计算事后资产配置的 Bias 统计量，以评估事后资产配置。

/*
获取因子的Bias时序统计量和获取个股的特质收益统计量
Input：
   retOut   getRetTable()函数返回的结果
   index_name  指数代码
   method   等权方法或者流通市值方法 'equal'，'float_market'
Output:
   Bias统计量
*/

事前资产配置指在实际收益数据可用之前，根据模型的预测和假设进行的资产配置。这个过程发生在投资决策之前，基于模型的预测结果和投资者的目标、约束条件等进行资产配置。根据已经由优化目标得到组合权重或是给定的组合权重，可以计算出指定组合的偏差统计量和 Q 统计量，观察指定资产配置组合权重的合理性或是评估优化权重的好坏。

我们基于 getPortfolioAccuracy 接口以评估事前资产配置组合。

/* getPortfolioAccuracy
计算资产组合的时序偏差统计量、Q-统计量
Input:                                      
         facTable          NULL   (Default)   (getFSLevelFactor函数返回的）全因子回归表
         retOut            NULL  (Default)    （getRetTable函数返回的）全因子收益表    
         st                2022.01.03(Default)                
         et                2023.01.02(Default)
         index_code        指定资产组合        '000300'、'399101'
         method            权重配比方法        "float_value"  流通市值加权，"equal" 等权
Output:  
         accuracyTable   
*/

假设以沪深 300 指数的等权资产组合为例，基于 getPortfolioAccuracy 函数计算事前预测资产组合的偏差，发现 B_n 统计量的均值为 0.335<1，但此处 Bias 结果与事前预测得到的组合权重有关，因此说明沪深 300 指数的等权组合对未来风险的可预测性较弱，建议重新配置资产组合。

/* 计算预测模型的bias、q统计量 */
index_code = '000300'
st = 2022.01.03 
et = 2023.01.02
outAccurary = getPortfolioAccuracy(st,et,facTable1,retOut,index_code,'equal')
outAccurary = outAccurary[2:]
baseline = take(1,(shape outAccurary)[0])
plot(table(outAccurary.bias_statistic,baseline),outAccurary.record_date, extras={multiYAxes: false})
mean(outAccurary) 

若要计算某个风格因子可通过get+风格因子名的函数获取。以 ABS 因子为例，可通过 getAbs() 函数计算得到，调用方法如下：

// Get raw style fatcors,some sample codes: 
getAbs(startTime = 2022.01.03,endTime = 2023.01.02)

// Get raw industry factors,some sample codes:
getIndustry(startTime = 2022.01.01,endTime = 2023.01.02,method = 'CITIC')
// Get Industry Factor Weights,some sample codes:
getIndustryWeighted(startTime = 2022.01.03,endTime = 2023.01.02,method = 'CITIC')
// Get weighted industry factors 
getIndustryFactor(startTime = 2022.01.03,endTime = 2023.01.02,method = 'CITIC')

RiskFactors 模块中的因子合成流程较为复杂，需要按照如下链路逐一处理

1. 基于 getAllFactors 获取所有因子
2. 基于 getGetRegTable 对原始因子表的缺失值进行处理
3. 基于 getFactorsValidation 对所有因子表进行有效性检验
4. 基于 getFSLevelFactor 合成因子

本文为便于用户进一步构建 RiskFactors 模型，用户可基于 prepareModelData 函数对上述流程统一调用。假设需要合成得到时间范围为 2022.01.03~2023.01.02、原始因子需要进行市值加权标准化、去极值处理、行业分类以中信一级 CITIC 行业分类为准、按照等权法合成的一级因子表，使用示例如下（若无需合成，指定 merge_level=NULL）：

st = 2022.01.03 
et = 2023.01.02
normlizing = true
scaling = true
decap = true
// also you can choose industry_method = 'SW_2021'
industry_method = 'CITIC'
industry_weighted = true
// Get ALL FIRST level of factors and factor Validatiion table
factorsValid,facTable1=prepareModelData(st,et,normlizing,scaling,decap,
                                        industry_method,industry_weighted,
                                        merge_method="equal",merge_level="F")

select * from facTable1 limit 100

x,y=plotFactorsValidation(factorsValid,"fsc")
plot(x,y,extras={multiYAxes: false},title = "因子fsc 月频时序图")

x,y=plotFactorsValidation(factorsValid,"ic")
plot(x,y,extras={multiYAxes: false},title = "因子ic月频时序图")

x,y=plotFactorsValidation(factorsValid,"t")
plot(x,y,extras={multiYAxes: false},title = "因子t_stat月频时序图")

基于 getRetTable 函数，用户可以基于因子表构建 RiskFactors 多因子模型，假设当市场收益率存在序列自相关、当市场关联密切且需要对个股特异性风险进行调整时，使用示例如下：

retOut1 = getRetTable(facTable1,adjust = true,shrink = true,eigenfactor = true)   

// the output of getRetTable function 
// the trait return covariance matrix for 12.30
retOut1.stock_risk[string(2022.12.30)] 
// risk factor covariance matrix for 12.30
retOut1.fac_risk[string(2022.12.30)] 
// R2
retOut1.R2 
// trait return
retOut1.res 
// t-stat
retOut1.tstat 
// factor returns
retOut1.fac_ret 
// bias statistic
retOut1.bias 

// R2 Monthly Frequency Timing Chart
plot(retOut1.R2.stR2,retOut1.R2.record_date,"𝑆𝑡𝑢𝑑𝑒𝑛𝑡𝑖𝑧𝑒𝑑 R2 月频时序图")

以最小风险为目标，假设预测得到的下一期个股收益率前20的股票为股票池，且无风险利率为0.05，则求得最优投资组合的使用示例如下：

// Initial screening stock1
optionCode = exec stock_code from getPredicOut(facTable1).predict_ret 
             order by return_day desc limit 20         

// Controlled return, minimized risk model
// Obtaining weight combinations
portWeight1 = getOptimizeWeights(facTable = facTable1,retOut = retOut1,st = st,et = et, 
                                  method ="minRiskControlRet",r = 0.05,
                                  optionCode = optionCode)  

若用户需要进行事后的因子投资组合评估，并绘制因子模型因子 Bias 统计量时序图、因子模型特异风险 Bias 统计量时序图，使用示例如下：

x,y=FactorCombinationAssessment(retOut = retOut1,plot_index="bias_stat")
plot(x,y,title = "因子模型因子Bias统计量时序图")

x,y = FactorCombinationAssessment(retOut = retOut1,plot_index="avg_bias_stat")
plot(x,y,extras={multiYAxes: false},title = "因子均值Bias统计量时序图")   

x,y = FactorCombinationAssessment(retOut = retOut1,plot_index="stock_bias_stat")
plot(x,y,extras={multiYAxes: false},title = "因子模型特异风险Bias统计量时序图") 

以沪深300指数为例，并以等权法构建指数，指数的简单资产配置评估使用示例如下：

x,y =AssetPortfolioAssessment(retOut = retOut1,index_code="000300",
                              index_weight="equal") 
plot(x,y,extras={multiYAxes: false})