初探多因子选股：基于Fama-Macbeth回归的因子分析框架 （附Python3代码）

Fama-Macbeth回归及因子统计

引言

本文介绍的因子统计方法基于1973年Fama和Macbeth为验证CAPM模型而提出的Fama-Macbeth回归，该模型现如今被广泛用被广泛用于计量经济学的panel data分析，而在金融领域在用于多因子模型的回归检验，用于估计各类模型中的因子暴露和因子收益（风险溢价）。

Fama-Macbeth与传统的截面回归类似，本质上也与是一个两阶段回归，不同的是它用了巧妙的方法解决了截面相关性的问题，从而得出更加无偏，相合的估计。

时间序列回归

Fama-Macbeth模型与传统截面回归相同，第一步都是做时间序列回归。在因子分析框架中，时间序列回归是为了获得个股在因子上的暴露。如果模型中的因子是 portfolio returns（即使用投资组合收益率作为因子，例如Fama-French三因子模型中的SMB，HML和市场因子），那么可以通过时间序列回归（time-series regression）来分析$E[R_i]$E[Ri​]和$\beta_i$βi​在截面上的关系。（本文举例的因子都是portfolio returns）

令$f_t$ft​为因子组合在t期的收益率，$R_{it}$Rit​为个股$i$i在t期的收益率，用$f_t$ft​对每只股票的$R_{it}$Rit​回归，即可得到每支股票的全样本因子暴露$\beta_i$βi​。
$R_{it}=\alpha_i+\beta_if_t+\varepsilon_{it},t=1,2,...,T \forall i$Rit​=αi​+βi​ft​+εit​,t=1,2,...,T∀i
也可滚动计算某个时间段的因子暴露$\beta_{it}$βit​,体现个股随市场的变化设置时间段长度为$period$period
$R_{ik}=\alpha_i+\beta_{it}f_k+\varepsilon_{ik},k=t-period,2,...,t \forall i$Rik​=αi​+βit​fk​+εik​,k=t−period,2,...,t∀i

截面回归

截面回归的第一步就是通过时间序列回归得到个股暴露，与Fama-Macbeth回归相同，第二步回归体现了传统截面回归和Fama-Macbeth回归最大的不同

对时序回归中回归式在时间序列上取均值，在$E[\varepsilon]=0$E[ε]=0的假设下可以得出：
$E[R_{i}]=\alpha_i+\beta_iE[f]$E[Ri​]=αi​+βi​E[f]
上式正是个股的期望收益与因子暴露在截面上的关系，截距$\alpha_i$αi​为个股的错误定价

那么便可通过截面回归找到因子的期望收益率$E[f]$E[f],方法是最小化个股定价错误$\alpha_i$αi​的平方和。对个股的的收益在时序上取均值得到个股期望收益$E[R_i]$E[Ri​]，用全样本的个股因子暴露对个股期望收益做无截距回归。
$E[R_{i}]=\beta_i\lambda+\alpha_i$E[Ri​]=βi​λ+αi​
回归残差$\alpha_i$αi​为个股的错误定价，$\lambda$λ为因子的期望收益率。

截面回归最大的缺陷在于忽略了截面上的残差相关性，使得OLS给出的标准误存在巨大的低估。

Fama-Macbeth回归

与截面回归相同，第一步都是通过时间序列回归得到因子暴露值，不同的是，第二步中，Fama-Macbeth在每个t上都做了一次无截距截面回归，：
$R_{it}=\beta_i\lambda_t+\alpha_{it},i=1,2,...,N \forall t$Rit​=βi​λt​+αit​,i=1,2,...,N∀t
上式中的$\beta_i$βi​为全样本$\beta$β，当然若使用滚动回归数据，也可以在不同截面的回归上使用对应时期的$\beta_{i,t}$βi,t​

Fama-Macbeth回归相当于在每个t上做一次独立的截面回归，这T次回归的参数取均值作为回归的估计值：
$\hat{\lambda}=\frac{1}{T} \sum_{t=1}^{T} \hat{\lambda}_{t},\hat{\alpha_i}=\frac{1}{T} \sum_{t=1}^{T} \hat{\alpha}_{it}$λ^=T1​t=1∑T​λ^t​,αi​^​=T1​t=1∑T​α^it​
上述方法的巧妙之处在于它把 T 期的回归结果当作 T 个独立的样本。参数的 standard errors 刻画的是样本统计量在不同样本间是如何变化的。在传统的截面回归中，我们只进行一次回归，得到$\lambda$λ和$\alpha_i$αi​的一个样本估计。而在 Fama-MacBeth 截面回归中，我们把T期样本点独立处理，得到 T 个$\lambda$λ和$\alpha_i$αi​的样本估计。

若使用全样本因子暴露$\beta_i$βi​进行估计，截面回归和Fama-Macbeth的估计结果相同，当使用滚动窗口进行估计时（Fama and MacBeth (1973)中作者使用了滚动窗口），截面回归和Fama-Macbeth回归会得到完全不同的估计结果。

Fama-Macbeth回归很好的解决了截面相关性的问题，但对于时间序列上的相关性仍然无力。

因子统计

Fama-Macbeth回归为本文所讲的因子统计框架提供了大量参数，包括每次截面回归的斜率$\lambda_t$λt​和每次回归系数的t值$tvalue_t$tvaluet​

图中的$beta_t$betat​就是每个截面上的$\lambda_t$λt​

Python3代码

本文将所有时序回归，截面回归和Fama-macbeth回归都封装在一个类里，方便调用。因为要进行多次的回归，最多N*T次，故没有使用第三方库，而是用OLS的矩阵解析式直接计算得到回归参数，测试出速度大概能比第三方库快3~4倍。代码已经尽笔者所能优化到了最快的速度，欢迎各位大佬搬运测试。

初始化

'''
@Time    : 2020/8/5 13:33
@Author  : hjz
@Email   : [e-mail]

'''
class Fama_regression:
    def __init__(self, return_data, factor_data, frequency='d'): 
     # return_data:T*N factor_data:T*1 frequency='d' , 'm' ,'y'
        self.T, self.N = return_data.shape[0], return_data.shape[1]
        self.stock = return_data.columns #股票池
        self.time = return_data.index   #时间
        self.return_data = return_data  # 股票收益率矩阵
        self.factor_data = factor_data  # 因子收益率序列
        if frequency=='d':              #频率（日，月，年）
            self.time_period=250
        elif frequency=='m':
            self.time_period=12
        else:
            self.time_period=1

时间序列回归（全样本）

    def time_series_regression_all_sample(self):  # 全样本时间序列回归
        def time_series_regression(Y):
            Y = Y.values
            Y[np.isnan(Y)] = 0
            X = self.factor_data.values
            constant = np.array([[1]] * len(self.time))
            X = np.hstack((constant, X))
            beta = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(Y).T[1]  # OLS矩阵求解式
            return beta
        self.factor_loading_allsample = self.return_data.apply(time_series_regression)
        return

时间序列回归（滚动窗口）

若参数Forecast为true表示预测，会将当期因子暴露记录在下一次，则之后的fama-macbeth回归中，当期因子暴露与下期收益回归，验证因子对收益的预测能力

    def time_series_regression_rolling(self, period=22,Forecast = False):# 滚动时间窗口时间序列回归，
        # period表示滚动窗口的长度
        # Forecast为true表示预测，用于验证因子对收益的预测能力，即后面的fama-macbeth回归用当期暴露对下期收益回归，则在此函数中将当期暴露记录在下期
        # Forecast为False表示后面的fama-macbeth回归用当期暴露对当期收益回归，表示对当期收益的归因
        if period > len(self.time):
            exit("Period is too large")
            return
        def time_series_regression(Y):
            if not hasattr(time_series_regression,'count'):
                time_series_regression.count=0
            time_series_regression.count+=1
            constant = np.array([[1]] * period)
            Y[np.isnan(Y)] = 0
            X = factor_data_temp.iloc[0:period, :].values
            X = np.hstack((constant, X))
            beta = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(Y).T[1]  # OLS矩阵求解式
            if time_series_regression.count==len(self.stock):
                factor_data_temp.drop(factor_data_temp.iloc[0, :].name, axis=0, inplace=True)
                time_series_regression.count = 0
            return beta
        factor_data_temp=self.factor_data
        rolling_data = self.return_data.rolling(window=period).apply(time_series_regression)
        if Forecast:
            rolling_data=rolling_data.iloc[period-1:-1]
            rolling_data.index = (self.time[period:])
        else:
            rolling_data=rolling_data.iloc[period-1:]
        self.df_factor_loading_rolling=rolling_data

        return

这一块的速度一直很慢，dataframe.rolling迭代效率太低，若各位大佬有更快的方法欢迎指教

截面回归

    def section_regression_without_intercept(self):  # 截面回归
        Er = self.return_data.apply(np.mean)
        N = len(Er)
        Y = np.array([Er.values.tolist()]).T
        X = np.array([self.factor_loading_allsample.values.tolist()]).T
        beta = (X.T.dot(Y)) / (X.T.dot(X))[0][0]
        epsilon = Y - beta * X
        T_test = beta / np.sqrt(epsilon.T.dot(epsilon) / (X.T.dot(X)) / (N - 1))[0][0]
        self.error_of_section_reg = pd.Series(epsilon.T[0], index=self.stock)
        self.beta_of_section_reg = beta
        self.tvalue_of_section_reg = T_test
        return

Fama-Macbeth回归

参数data_type表示使用全样本因子暴露还是滚动窗口因子暴露

    def fama_macbeth_regression_without_intercept(self,data_type='rolling'):  # Fama-macbeth回归
        #data_type为rolling（滚动回归数据），allsample（全样本回归数据）
        def section_regression_epsilion(Y):
            if data_type=='rolling':
                X = self.df_factor_loading_rolling.loc[Y.name]
            else:
                X = self.factor_loading_allsample
            N = len(X)
            X = np.array([X.values.tolist()]).T
            Y = np.array([Y.values.tolist()]).T
            beta = ((X.T.dot(Y)) / (X.T.dot(X)))[0][0]
            epsilon = Y - beta * X
            epsilon = epsilon.T[0]
            return pd.Series(epsilon, index=self.stock)

        def section_regression_beta(Y):
            if data=='rolling':
                X = self.df_factor_loading_rolling.loc[Y.name]
            else:
                X = self.factor_loading_allsample
            N = len(X)
            X = np.array([X.values.tolist()]).T
            Y = np.array([Y.values.tolist()]).T
            beta = ((X.T.dot(Y)) / (X.T.dot(X)))[0][0]
            epsilon = Y - beta * X
            T_test = beta / np.sqrt(epsilon.T.dot(epsilon) / (X.T.dot(X)) / (N - 1))[0][0]
            return pd.Series([beta, T_test], index=['beta', 'tvalue'])

        time_list = self.df_factor_loading_rolling.index
        return_data = self.return_data.loc[time_list]
        self.epsilon_mat = return_data.apply(section_regression_epsilion, axis=1)
        self.epsilon_mean = self.epsilon_mat.apply(np.mean)
        self.beta_tvalue = return_data.apply(section_regression_beta, axis=1)
        self.beta_fama, self.tvalue_fama = self.beta_tvalue['beta'], self.beta_tvalue['tvalue']
        return

因子计算及显示

 def compute_factor_characteristic(self):
        def autocor(X):
            if not hasattr(autocor,'last'):
                autocor.last=X.values.copy()
                return 0
            else:
                result=np.corrcoef(X.values,autocor.last)
                autocor.last=X.values.copy()
                return result[0][1]
        def IC_rank(X):
            R=self.return_data.loc[X.name]
            mat=pd.DataFrame([X,R]).T
            return mat.corr('spearman').values[0][1]

        self.autocor = self.df_factor_loading_rolling.apply(autocor, axis=1)[1:]
        self.IC_rank=self.df_factor_loading_rolling.apply(IC_rank,axis=1)
        self.IC_rank_mean=self.IC_rank.mean()
        self.IC_IR_rank=self.IC_rank.std()
        self.factor_return_annual=self.beta_fama.mean()*self.time_period
        self.factor_vol_annual = self.beta_fama.std() * np.sqrt(self.time_period)
        self.factor_sharpe_ratio=self.factor_return_annual/self.factor_vol_annual
        self.factor_tvalue=self.beta_fama.mean()/(self.beta_fama.mean()*np.sqrt(len(self.beta_fama)))
        self.mean_tvalue=self.tvalue_fama.mean()
        self.mean_abs_tvalue=np.mean(np.abs(self.tvalue_fama))
        self.tvalue_morethan2=(self.tvalue_fama.abs()>2).sum()/len(self.tvalue_fama)

    def show_factor_characteristic(self):
        print("因子截面相关性: ",round(self.autocor.mean(),2))
        print("因子IC：",round(self.IC_rank_mean,2))
        print("因子IC_IR：", round(self.IC_IR_rank, 2))
        print("因子年化收益：", round(self.factor_return_annual, 2))
        print("因子年化波动率：", round(self.factor_vol_annual, 2))
        print("因子夏普比率: ",round(self.factor_sharpe_ratio, 2))
        print("因子t值: ", round(self.factor_tvalue, 2))
        print("平均t值: ", round(self.mean_tvalue, 2))
        print("平均绝对t值: ", round(self.mean_abs_tvalue, 2))
        print("绝对t值>2占比: ", round(self.tvalue_morethan2, 2))
        return

测试结果

本文选用沪深300成分股2020年的日数据对市场因子beta进行测试
输出结果：

因子截面相关性:  -0.04
因子IC： -0.0
因子IC_IR： 0.06
因子年化收益： 0.46
因子年化波动率： 0.25
因子夏普比率:  1.82
因子t值:  0.09
平均t值:  1.23
平均绝对t值:  9.93
绝对t值>2占比:  0.83

时间开销（单位秒）：

时间开销（单位：秒）：
时序回归全样本：  0.04852179628497339
时序回归全滚动窗口：  11.162792468957404
截面回归：  0.039669358541182476
Fama-Macbeth回归：  0.1491620773626554
因子统计计算：  1.4921370042123563

在滚动窗口时序回归上要花费大量时间，我测试过很多方法去替代dataframe.rolling，但结果都不理想，欢迎各位大佬提出改进检验。

参考

股票多因子模型的回归检验——石川
《长江证券-金融工程专题-覃川桃郑起-高频因子（三）：高频因子研究框架》