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

方法（median absolute deviation 方法），该方法使用中位数，较于传统的均值标准差具有稳健（robust）性质，不会被极端异常值影响结果，是对减均值后标准差处理的改进方法。 方法相对于其他离群值检测方法的优点在于它对离群值的鲁棒性较强。本文中 方法对应接口函数 winsorized 。

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

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

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

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

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

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

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

本文市值加权标准化对应接口函数 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，本文选择采用均值填充以处理缺失值。

/* 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，加权最小二乘）以为因变量，以风格因子 s 的因子载荷为自变量，预测（风格因子的收益率）。

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

  此处使用 以为因变量，以行业因子 i 的因子载荷为自变量，预测（行业因子的收益率）

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

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

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

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

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

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

Information Coefficient（IC）值是因子对下期收益率的预测与下期实际收益率的相关性，IC 代表的是预测值和实现值之间的相关性，通常用于评价预测能力（即选股能力）。在实际计算中，因子 𝑘 的 IC 值一般是指个股第 𝑇 期在因子 𝑘 上的暴露度与 𝑇 +1 期的收益率的相关系数。因子 IC 值反映的是个股下期收益率和本期因子暴露度的线性相关程度，表示使用该因子进行收益率预测的稳健性。IC 的计算方式有两种： IC、rank IC。本文在 getFactorsValidation 函数基于 Spearman 相关系数计算 Rank IC 值。

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

本文的单因子有效性检验的对应接口函数 getFactorsValidation。该函数通过输入计算得到的全因子表，给出了每一个单因子模型的单因子收益，t 统计量，拟合优度， 指标和 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)

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

  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 月频时序图")

  
  
  

  大部分处于 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）是一个常见的将先验和样本估计值结合起来的手段。进行贝叶斯收缩调整的主要原因是，基于 Barram 模型使用样本内数据计算出的特异性波动率在样本外的持续性很差，极高或极低的个股波动率可能会出现均值回归的情况，那么波动率会被高估或低估。

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

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

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

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

  其中 Barra 把所有个股按照市值分成十档，表示对应的市值档位，为同处在市值档位下的 n 支股票的特异性收益率的波动，表示 n 支股票按照市值计算的权重。

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

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

  
  
  

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

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

• 行业中性

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

  
  
  
• 风格因子中性

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

  
  
  
• 现金中性

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

  
  
  
• 最小收益预测

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

  
  
  

本文给出行业基准中性和风格基准中性下的控制最小预测收益并最小组合风险方法，对应接口函数 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 full original fatcors    查看原始风格三级因子
getAbs(startTime = 2022.01.03,endTime = 2023.01.02)
getAcf(method='TTM',startTime = 2022.01.03,endTime = 2023.01.02)                
//'TTM'(Default)
getAcf(method='LYR',startTime = 2022.01.03,endTime = 2023.01.02)
getVsal(method='TTM',startTime = 2022.01.03,windows = 1825,endTime = 2023.01.02)
getVsal(method='LYR',startTime = 2022.01.03,windows = 1825,endTime = 2023.01.02)   
//'LYR'(Default)
getVern(method='TTM',startTime = 2022.01.03,windows = 1825,endTime = 2023.01.02)    
getVern(method='LYR',startTime = 2022.01.03,windows = 1825,endTime = 2023.01.02)   
//'LYR'(Default)
getVflo(method='TTM',startTime = 2022.01.03,windows = 1825,endTime = 2023.01.02)
getVflo(method='LYR',startTime = 2022.01.03,windows = 1825,endTime = 2023.01.02)  
 //'LYR'(Default)
getCetop(method='TTM',startTime = 2022.01.03,endTime = 2023.01.02)              
//'TTM'(Default)
getCetop(method='LYR',startTime = 2022.01.03,endTime = 2023.01.02)              
getEtop(method='TTM',startTime = 2022.01.03,endTime = 2023.01.02)               
//'TTM'(Default)
getEtop(method='LYR',startTime = 2022.01.03,endTime = 2023.01.02)
getEgro(method='TTM',startTime = 2022.01.03,windows = 1825,endTime = 2023.01.02)
getEgro(method='LYR',startTime = 2022.01.03,windows = 1825,endTime = 2023.01.02)   
//'LYR'(Default)
getSgro(method='TTM',startTime = 2022.01.03,windows = 1825,endTime = 2023.01.02)
getSgro(method='LYR',startTime = 2022.01.03,windows = 1825,endTime = 2023.01.02)   
//'LYR'(Default)
getAgro(startTime = 2022.01.03,windows = 1825,endTime = 2023.01.02)
getIgro(startTime = 2022.01.03,windows = 1825,endTime = 2023.01.02)
getMlev(startTime = 2022.01.03,windows = 365,endTime = 2023.01.02)
getDtoa(startTime = 2022.01.03,windows = 365,endTime = 2023.01.02)
getBlev(startTime = 2022.01.03,windows = 365,endTime = 2023.01.02)   
getStom(startTime = 2022.01.03,windows = 21,endTime = 2023.01.02)
getStoq(startTime = 2022.01.03,windows = 63,endTime = 2023.01.02)
getStoa(startTime = 2022.01.03,windows = 365,endTime = 2023.01.02)
getAtvr(startTime = 2022.01.03,windows = 252,halfLife = 63,endTime = 2023.01.02)
getLtrstr(startTime = 2022.01.03,windowsReturn = 1040,windowsAvg = 11,lagged = 273, halfLife = 260,endTime = 2023.01.02)
getLthalpha(startTime = 2022.01.03,windowsReturn = 1040,windowsAvg = 11, lagged = 273,halfLife = 260,endTime = 2023.01.02)
getMidcap(startTime = 2022.01.03,endTime = 2023.01.02)
getRstr(startTime = 2022.01.03,windowsReturn = 252,windowsAvg = 11,lagged = 11, halfLife = 126,endTime = 2023.01.02)
getHalpha(startTime = 2022.01.03,windowsReturn = 252,windowsAvg = 11,lagged = 11, halfLife = 126,endTime = 2023.01.02)
getAto(method='TTM',startTime = 2022.01.03,endTime = 2023.01.02)                
//'TTM'(Default)
getAto(method='LYR',startTime = 2022.01.03,endTime = 2023.01.02)
getRoa(method='TTM',startTime = 2022.01.03,endTime = 2023.01.02)               
 //'TTM'(Default)
getRoa(method='LYR',startTime = 2022.01.03,endTime = 2023.01.02)
getBtop(startTime = 2022.01.03,endTime = 2023.01.02)
getDtop(startTime = 2022.01.03,windows = 365,endTime = 2023.01.02)
getHbeta(startTime = 2022.01.03,windowsReturn = 252,halfLife = 63, endTime = 2023.01.02)
getHsigma(startTime = 2022.01.03,windowsReturn = 252,halfLife = 63, endTime = 2023.01.02)
getDastd(startTime = 2022.01.03,windowsReturn = 252,halfLife = 42, endTime = 2023.01.02)
getCmra(startTime = 2022.01.03,windowsReturn = 365,endTime = 2023.01.02)
getLncap(startTime = 2022.01.03,endTime = 2023.01.02)
getCxgro(startTime = 2022.01.03,windows = 5,endTime = 2023.01.02)
getGp(method='TTM',startTime = 2022.01.03,windows = 365,endTime = 2023.01.02)
getGp(method='LYR',startTime = 2022.01.03,windows = 365,endTime = 2023.01.02)    
 //'LYR'(Default)
getGpm(method='TTM',startTime = 2022.01.03,windows = 365,endTime = 2023.01.02)
getGpm(method='LYR',startTime = 2022.01.03,windows = 365,endTime = 2023.01.02)    
//'LYR'(Default)
getEm(startTime = 2022.01.03,windows = 365,endTime = 2023.01.02)

有了计算单个风格因子的函数后，调用 getIndustry 等函数获得行业因子：

// 获取原始行业因子
getIndustry(startTime = 2022.01.01,endTime = 2023.01.02,method = 'SW_2021') 
getIndustry(startTime = 2022.01.01,endTime = 2023.01.02,method = 'CITIC')
// 获取行业因子权重
getIndustryWeighted(startTime = 2022.01.03,endTime = 2023.01.02,method = 'SW_2021')
getIndustryWeighted(startTime = 2022.01.03,endTime = 2023.01.02,method = 'CITIC')
// 获取加权后的行业因子
getIndustryFactor(startTime = 2022.01.03,endTime = 2023.01.02,method = 'SW_2021')
getIndustryFactor(startTime = 2022.01.03,endTime = 2023.01.02,method = 'CITIC')

在计算得到所有风格因子后，可以调用 standardized、winsorized 对风格因子进行标准化、去极值处理。

//factors use for factorValidation 去极值、标准化、市场中性   
//在同种类别因子中，如若包含ttm或者lyr，只能选取其中一种纳入模型 
abs = standardized(winsorized(getAbs(startTime = 2022.01.03,endTime = 2023.01.02)),
                   adjusted = true)
acf_ttm = standardized(winsorized(getAcf(method='TTM',startTime = 2022.01.03, endTime = 2023.01.02)),adjusted = true)   //'TTM'(Default)
acf_lyr = standardized(winsorized(getAcf(method='LYR',startTime = 2022.01.03, endTime = 2023.01.02)),adjusted = true)
vsal_ttm = standardized(winsorized(getVsal(method='TTM',startTime = 2022.01.03, windows = 1825,endTime = 2023.01.02)),adjusted = true)
vsal_lyr = standardized(winsorized(getVsal(method='LYR',startTime = 2022.01.03, windows = 1825,endTime = 2023.01.02)),adjusted = true) 
                   //'LYR'(Default)
vern_ttm = standardized(winsorized(getVern(method='TTM',startTime = 2022.01.03, windows = 1825,endTime = 2023.01.02)),adjusted = true)    
vern_lyr = standardized(winsorized(getVern(method='LYR',startTime = 2022.01.03, windows = 1825,endTime = 2023.01.02)),adjusted = true)   
                  //'LYR'(Default)
vflo_ttm = standardized(winsorized(getVflo(method='TTM',startTime = 2022.01.03, windows = 1825,endTime = 2023.01.02)),adjusted = true)
vflo_lyr = standardized(winsorized(getVflo(method='LYR',startTime = 2022.01.03, windows = 1825,endTime = 2023.01.02)),adjusted = true)
                    //'LYR'(Default)
cetop_ttm = standardized(winsorized(getCetop(method='TTM',startTime = 2022.01.03, endTime = 2023.01.02)),adjusted = true)             
                   //'TTM'(Default)
cetop_lyr = standardized(winsorized(getCetop(method='LYR',startTime = 2022.01.03, endTime = 2023.01.02)),adjusted = true)             
etop_ttm = standardized(winsorized(getEtop(method='TTM',startTime = 2022.01.03, endTime = 2023.01.02)),adjusted = true)              
                    //'TTM'(Default)
etop_lyr = standardized(winsorized(getEtop(method='LYR',startTime = 2022.01.03, endTime = 2023.01.02)),adjusted = true)
egro_ttm = standardized(winsorized(getEgro(method='TTM',startTime = 2022.01.03, windows = 1825,endTime = 2023.01.02)),adjusted = true)
egro_lyr = standardized(winsorized(getEgro(method='LYR',startTime = 2022.01.03, windows = 1825,endTime = 2023.01.02)),adjusted = true)
                    //'LYR'(Default)
sgro_ttm = standardized(winsorized(getSgro(method='TTM',startTime = 2022.01.03, windows = 1825,endTime = 2023.01.02)),adjusted = true)
sgro_lyr = standardized(winsorized(getSgro(method='LYR',startTime = 2022.01.03, windows = 1825,endTime = 2023.01.02)),adjusted = true)  
                   //'LYR'(Default)
agro = standardized(winsorized(getAgro(startTime = 2022.01.03,windows = 1825, endTime = 2023.01.02)),adjusted = true)
igro = standardized(winsorized(getIgro(startTime = 2022.01.03,windows = 1825, endTime = 2023.01.02)),adjusted = true)
mlev = standardized(winsorized(getMlev(startTime = 2022.01.03,windows = 365, endTime = 2023.01.02)),adjusted = true)
dtoa = standardized(winsorized(getDtoa(startTime = 2022.01.03,windows = 365, endTime = 2023.01.02)),adjusted = true)
blev = standardized(winsorized(getBlev(startTime = 2022.01.03,windows = 365, endTime = 2023.01.02)),adjusted = true)
stom = standardized(winsorized(getStom(startTime = 2022.01.03,windows = 21, endTime = 2023.01.02)),adjusted = true)
stoq = standardized(winsorized(getStoq(startTime = 2022.01.03,windows = 63, endTime = 2023.01.02)),adjusted = true)
stoa = standardized(winsorized(getStoa(startTime = 2022.01.03,windows = 365, endTime = 2023.01.02)),adjusted = true)
atvr = standardized(winsorized(getAtvr(startTime = 2022.01.03,windows = 252, halfLife = 63,endTime = 2023.01.02)),adjusted = true)
ltrstr = standardized(winsorized(getLtrstr(startTime = 2022.01.03, windowsReturn = 1040,windowsAvg = 11,lagged = 273, halfLife = 260,endTime = 2023.01.02)),adjusted = true)
lthalpha = standardized(winsorized(getLthalpha(startTime = 2022.01.03, windowsReturn = 1040,windowsAvg = 11,lagged = 273, halfLife = 260,endTime = 2023.01.02)),adjusted = true)
midcap = standardized(winsorized(getMidcap(startTime = 2022.01.03, endTime = 2023.01.02)),adjusted = true)
rstr = standardized(winsorized(getRstr(startTime = 2022.01.03, windowsReturn = 252,windowsAvg = 11,lagged = 11, halfLife = 126,endTime = 2023.01.02)),adjusted = true)
halpha = standardized(winsorized(getHalpha(startTime = 2022.01.03, windowsReturn = 252,windowsAvg = 11,lagged = 11, halfLife = 126,endTime = 2023.01.02)),adjusted = true)
ato_ttm = standardized(winsorized(getAto(method='TTM',startTime = 2022.01.03, endTime = 2023.01.02)),adjusted = true)   //'TTM'(Default)
ato_lyr = standardized(winsorized(getAto(method='LYR',startTime = 2022.01.03, endTime = 2023.01.02)),adjusted = true)  
roa_ttm = standardized(winsorized(getRoa(method='TTM',startTime = 2022.01.03, endTime = 2023.01.02)),adjusted = true)   //'TTM'(Default)
roa_lyr = standardized(winsorized(getRoa(method='LYR',startTime = 2022.01.03, endTime = 2023.01.02)),adjusted = true)
btop = standardized(winsorized(getBtop(startTime = 2022.01.03, endTime = 2023.01.02)),adjusted = true)
dtop = standardized(winsorized(getDtop(startTime = 2022.01.03,windows = 365, endTime = 2023.01.02)),adjusted = true)
hbeta = standardized(winsorized(getHbeta(startTime = 2022.01.03, windowsReturn = 252,halfLife = 63, endTime = 2023.01.02)),adjusted = true)
hsigma = standardized(winsorized(getHsigma(startTime = 2022.01.03, windowsReturn = 252, halfLife = 63,endTime = 2023.01.02)), adjusted = true)
dastd = standardized(winsorized(getDastd(startTime = 2022.01.03, windowsReturn = 252,halfLife = 42, endTime = 2023.01.02)),adjusted = true)
cmra = standardized(winsorized(getCmra(startTime = 2022.01.03,windowsReturn = 365, endTime = 2023.01.02)),adjusted = true)
lncap = standardized(winsorized(getLncap(startTime = 2021.12.01, endTime = 2023.01.02)),adjusted = false)
cxgro = standardized(winsorized(getCxgro(startTime = 2022.01.03, windows = 1825, endTime = 2023.01.02)),adjusted = true)
gp_ttm = standardized(winsorized(getGp(method='TTM',startTime = 2022.01.03, windows = 365,endTime = 2023.01.02)),adjusted = true)
gp_lyr = standardized(winsorized(getGp(method='LYR',startTime = 2022.01.03, windows = 365,endTime = 2023.01.02)),adjusted = true) 
                   //'LYR'(Default)
gpm_ttm = standardized(winsorized(getGpm(method='TTM',startTime = 2022.01.03, windows = 365,endTime = 2023.01.02)),adjusted = true)
gpm_lyr = standardized(winsorized(getGpm(method='LYR',startTime = 2022.01.03, windows = 365,endTime = 2023.01.02)),adjusted = true)  
                   //'LYR'(Default)
em = standardized(winsorized(getEm(startTime = 2022.01.03,windows = 365, endTime = 2023.01.02)),adjusted = true)

接下来，批量通过调用函数并行计算所有因子，包括风格因子和行业因子（如果需要查看单个风格因子或行业因子，推荐使用上述的单个因子计算函数）：

/* Barra 因子+模型全流程 */
st = 2012.12.01
et = 2023.01.02
normlizing = true
scaling = true
decap = false
industry_method = 'CITIC'
industry_weighted = false


// 获取所有因子、包括行业因子、风格因子
Factors = getAllFactors(st=st,et =et, normlizing = true,scaling = true,
                   decap = false,industry_method = 'CITIC',
                   industry_weighted = false)
Factors = getAllFactors(st=st,et =et, normlizing = true,scaling = true,
                   decap = false,industry_method = 'SW_2021',
                   industry_weighted = false)
select * from Factors limit 100

接下来对三级因子进行单因子检验，并画出相应的检验指标：（具体画图结果可参考本教程 2.4.5 的结果）

st = 2022.01.03 
et = 2023.01.02
normlizing = true
scaling = true
decap = true
industry_method = 'CITIC'
industry_weighted = true
Factors = getAllFactors(st= st,et = et, normlizing = normlizing,
                        scaling = scaling,decap = decap,
                        industry_method = industry_method,
                        industry_weighted = industry_weighted)

select * from Factors limit 100

// 对原始因子宽表的缺失值进行处理
fTable = getRegTable(factorsTable = true,tbName = Factors,st= st,et = et,
                    normlizing = normlizing ,scaling = scaling ,decap = decap,
                    industry_method = industry_method,
                    industry_weighted = industry_weighted)

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

update factorsValid set tstat = abs(tstat)

// 计算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 月频时序图")

// 计算ic
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 月频时序图")

// 计算tstat
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 月频时序图")

利用已经计算好的三级因子合成对应的一级、二级因子：

// 合成一级、二级因子 
//  一级因子和二级因子的json，可根据检验结果修改关系,若有需求，可自行添加其他因子
firstFactors = {
        "Quality":["Earnings_Quality","Earnings_Variability","
        Investment_Quality","Leverage","Profitability"],
        "Value":["Btop","Earning_Yield","Long_Term_Reversal"],
        "Growth":"Growth",
        "Liquidity":"Liquidity",
        "Volatility":["Beta","Residual_Volatility"],
        "Size":["Size","Mid_Cap"],
        "Momentum":"Momentum",
        "Dividend_Yield":"Dividend_Yield"
}   
    
// 二级因子和三级因子的json，可根据检验结果修改关系，三级因子一般需要选择ttm和lyr种类的因子
// 三级因子字典表
/* ['abs','acf_ttm','acf_lyr','vsal_ttm','vsal_lyr','vern_ttm','vern_lyr',
   'vflo_ttm','vflo_lyr','cetop_ttm','cetop_lyr','etop_ttm','etop_lyr',
   'egro_ttm','egro_lyr','sgro_ttm','sgro_lyr','agro','igro','mlev','dtoa',
   'blev','stom','stoq','stoa','atvr','ltrstr','lthalpha','midcap',
   'rstr','halpha','ato_ttm','ato_lyr','roa_ttm','roa_lyr','btop',
   'dtop','hbeta','hsigma','dastd','cmra','lncap','cxgro','gp_ttm',
   'gp_lyr','gpm_ttm','gpm_lyr','em'] */
secondFactors = {
    "Earnings_Quality":["abs","acf_ttm"],
    "Earnings_Variability":["vsal_ttm","vern_ttm","vflo_ttm"],
    "Investment_Quality":["agro","igro","cxgro"],
    "Leverage":["mlev","dtoa","blev"],
    "Profitability":["ato_ttm","gp_ttm","gpm_ttm","roa_ttm"],
    "Btop":"btop",
    "Earning_Yield":["cetop_ttm","etop_ttm","em"],
    "Long_Term_Reversal":["ltrstr","lthalpha"],
    "Growth":["egro_ttm","sgro_ttm"],
    "Liquidity":["stom","stoq","stoa","atvr"],
    "Beta":"hbeta",
    "Residual_Volatility":["hsigma","dastd","cmra"],
    "Size":"lncap",
    "Mid_Cap":"midcap",
    "Momentum":["rstr","halpha"],
    "Dividend_Yield":"dtop"
}

select count(*) from fTable
// 针对原始因子表填补缺失值
fTable = getRegTable(factorsTable = true,tbName = Factors,
                     st= st,et = et, normlizing = normlizing ,
                    scaling = scaling ,decap = decap,
                    industry_method = industry_method,
                    industry_weighted = industry_weighted)
// 原始因子有效性检验(factor_return、tstat、R2、aR2、stR2、fsc、IC)
factorsValid = getFactorsValidation(factorsTable = true,tbName = Factors,
                                    st = st,et = et , normlizing = normlizing,
                                    scaling = scaling,decap = decap,
                                    industry_method = industry_method,
                                    industry_weighted = industry_weighted)

select avg(R2)  from factorsValid group by valueType
select avg(aR2)  from factorsValid group by valueType
select avg(abs(tstat))  from factorsValid group by valueType

// 得到最终的一级、二级、因子（目前支持等权法、信息比率法、信息系数信息比率法）
facTable1 = getFSLevelFactor(fTable,factorsValid,firstFactors,secondFactors,false, 
                             "equal",level = "F")
facTable1 = getFSLevelFactor(fTable,factorsValid,firstFactors,secondFactors,false, 
                             "ir",level = "F")
facTable1 = getFSLevelFactor(fTable,factorsValid,firstFactors,secondFactors,false,
                              "ic_ir",level = "F")
facTable2 = getFSLevelFactor(fTable,factorsValid,firstFactors,secondFactors,true, 
                             "equal",level = "S")
facTable2 = getFSLevelFactor(fTable,factorsValid,firstFactors,secondFactors,true, 
                             "ir",level = "S")
facTable2 = getFSLevelFactor(fTable,factorsValid,firstFactors,secondFactors,true, 
                             "ic_ir",level = "S")

facTable1.columnNames()
select * from facTable1 limit 100

计算合成后的因子收益率和实现预测收益率：

// 一级因子收益率回归  
retOut1 = getRetTable(facTable1,adjust = true)                                    
// adjust采用Newey-West协方差调整、当市场收益率存在序列自相关时采用此方法
retOut1 = getRetTable(facTable1,adjust = false,shrink = false)                    
// shrink采用贝叶斯收缩调整特异性风险、推荐使用
retOut1 = getRetTable(facTable1,adjust = true,shrink = true)                      
// 特异性风险采用bayesian shrinkage
retOut1 = getRetTable(facTable1,adjust = true,shrink = true,eigenfactor = true)   
// 因子风险采用eigenfactor adjust,当市场关联密切时采用此方法
retOut1 = getRetTable(facTable1,adjust = true,shrink = true,eigenfactor = true)   
// 综上，推荐使用
retOut1 = getRetTable(facTable1,adjust = false,shrink = true,eigenfactor = true)  
// 综上，推荐使用
retOut1 = getRetTable(facTable1,adjust = false,shrink = true,eigenfactor = true)  
// 综上，推荐使用
retOut1 = getRetTable(fTable,adjust = false,shrink = true,eigenfactor = false)    
// 综上，推荐使用
undef(`retOut)
retOut = getRetTable(facTable1,adjust = true,shrink = false ,eigenfactor = false) 
 // 综上，推荐使用


// 例如
retOut1.stock_risk[string(2022.12.30)]   // 12.30的特质收益协方差矩阵
retOut1.fac_risk[string(2022.12.30)]     // 12.30的风险因子协方差矩阵
retOut1.R2                               // R2
retOut1.res                              // 特质收益
retOut1.tstat                            // t-stat
retOut1.fac_ret                          // 因子收益
retOut1.bias                             // bias统计量

select `All_Factors_R2,mean(stR2) from retOut.R2  where record_date 
between 2022.01.01 and 2023.01.01

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

retOut2 = getRetTable(facTable2)     // 二级因子收益率回归,相应可以得到对应模型检验结果
retOut3 = getRetTable(fTable)        // 三级因子收益率回归,相应可以得到对应模型检验结果

predictOut = getPredicOut(facTable1)         
// 结果较多，如需查看，推荐持久化后再查看结果，或者可以取出来部分较少的结果，除去预测收益外
pr = select * from predictOut.predict_ret    // 利用本期因子暴露预测最后一期的收益率
predictOut.R2                               // 预测模型R2
predictOut.res                              // 预测模型特质收益
predictOut.tstat                            // 预测模型t-stat
predictOut.fac_ret                          // 预测模型因子收益
predictOut.bias                             // 预测模型bias统计量

利用预测收益控制最小预测收益率，最小化风险，求解优化目标（需要注意，如果是模拟数据的话，由于模拟数据不满足一些模型条件，该步骤之后调用函数可能报错）：

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 统计量，需要给出月频的资产组合优化权重模型的结果。如下例子，模型不变，例如都是控制收益最小化风险，利用前 n-1 期数据预测第 n 期，给出第 n 期的统计量，如果需要获取每一期的统计量、每一期对模型给出估计后，都应该把前 t 期的 tbName、 facTable、 retOut结果代入到该函数（设为第 t 期），得到第 t 期的统计量结果。

/* 计算预测模型的bias、q统计量 */
// 真实数据
index_code = '000300'
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)

事后资产配置评估，计算 Barra 模型的 bias、q 统计量：

/*资产配置评估
*/
// 获取所有因子的时序bias统计量的值和所有个股的时序bias统计量
biasOut = getFacSpeacialBias(retOut)
// 因子bias
tmpfBias = select bias_stat from biasOut.fac_bias pivot by record_date,valueType
tbfBias = sql(select = sqlCol(tmpfBias.columnNames()[1:9]),from = tmpfBias).eval()

plot(tbfBias,tmpfBias.record_date,extras={multiYAxes: false},
      title = "因子模型因子Bias统计量时序图")
plot(tbfBias,tmpfBias.record_date,extras={multiYAxes: false})
code0 = parseExpr("rowAvg("+ concat(tmpfBias.columnNames()[1:],',') + ")")
avgfBias =  sql(select = sqlColAlias(code0,'avg_bias_stat'),from = tmpfBias).eval()
plot(avgfBias,tmpfBias.record_date,extras={multiYAxes: false},
      title = "因子均值Bias统计量时序图")
plot(avgfBias,tmpfBias.record_date,extras={multiYAxes: false})

// 个股特异bias
tmpsBias = select mean(bias_stat) from biasOut.stock_bias group by record_date
plot(tmpsBias.avg_bias_stat,tmpsBias.record_date,extras={multiYAxes: false},
      title = "因子模型特异风险Bias统计量时序图")

/* 简单资产配置评估 */
// 计算指数配置的Bias统计量
tmpIndexbiasbn = getFacSpeacialBias(retOut,'000300','equal').stock_bias
tmpIndexBias = select wavg(bias_stat,weight) from tmpIndexbiasbn group by record_date
plot(tmpIndexBias.wavg_bias_stat,tmpIndexBias.record_date,extras={multiYAxes: false})

tmpIndexbiasbn = getFacSpeacialBias(retOut,'000300','float_market').stock_bias
tmpIndexBias = select wavg(bias_stat,weight) from tmpIndexbiasbn group by record_date

plot(tmpIndexBias.wavg_bias_stat,tmpIndexBias.record_date,extras={multiYAxes: false})