Fama三因子模型应用及python代码

Fama-French三因子模型原理

自从资本资产定价模型（CAPM）问世以后，许多学者对其进行了实证研究，如Black和Scholes（1972）及Fama(1973)的检验证明，对1969年以前的数据而言，资本资产定价模型是有效的，而对之后的股市数据，却缺乏说服力。在横截面数据里，股票的平均收益和市场β相关性很低，因而更多影响股票收益的因素亟待发掘。

Fama和French(1993)开始对美国股市中股票收益率的决定因素进行了全面性的研究分析，他们发现单独使用Beta值或者分别与市值、P/E比、杠杆比例、B/M比结合在一起来研究股票收益率时，Beta的解释能力很弱。而市值、P/E比、杠杆比例、B/M比各因子单独用来解释收益率时，每个因子的解释能力都很强，当把这些因子组合起来时，市值、B/M比会弱化P/E比、杠杆比例的解释能力。

因此Fama和French从可以解释股票收益率的众多因素中提取出了三个重要的影响因子，即市场风险溢酬因子、市值因子和账面市值比因子B/M Ratio，仿照CAPM模型用这三个因子建立起来一个线性模型来解释股票的收益率，提出了著名的三因子模型（Fama and French Three Factor Model，简称FF3）。

三因子模型中的3个因子均为投资组合的收益率：市场风险溢酬因子对应了（市场投资组合的收益率减去无风险利率），市值因子对应了做多市值较小的公司与做空市值较大的公司的投资组合带来的收益率，账面市值比因子对应的是做多高BM公司、做空低BM公司的投资组合带来的收益率。三因子模型的形式为：

Rit-rf=αi+β1iRmt-rf+β2iSMBt+β3iHMLt+εit

其中Rit代表资产收益率，rf代表无风险收益率；Rit-rf为超额市场收益率；SMBt代表市值规模因子，HMLt代表账面市值比因子，β1i、β2i、β3i分别为Rit-rf、SMBt、HMLt的系数，εit为残差项，αi为截距项。

理论假设

Fama-French三因子模型，是以“有限理性”理论假设有基础，并在此基础上提出若干基本假定：

(1)存在着大量投资者；

(2)所有投资者都在同一证券持有期计划自己的投资资产组合；

(3)投资者投资范围仅限于公开金融市场上交易的资产；

(4)不存在证券交易费用(佣金和服务费用等)及税赋；

(5)投资者们对于证券回报率的均值、方差及协方差具有相同的期望值；

(6)所有投资者对证券的评价和经济局势的看法都一致。

三因子的构建方法

FF三因子模型提出，根据上市公司的市值，按照其大小值进行排序并分为两组，记为S、B，S、B分别表示为小市场规模股和大市场规模股。然后再根据年末上市公司账面市值比，按照33%、33%、33%的比例排序，记为L、M、H，L、M、H分别为低价值、中等价值、高价值；最后即可得到股票交叉组合，并通过加权平均（以总市值为权重）计算它们的月收益率{SL,SM,SH,BL,BM,BH}。通过以上6个组合的月收益率数据即可构造出市值规模因子（SMB），具体计算公式如下：

SMBt=(SLt+SMt+SHt)3-(BLt+BMt+BHt)3

通过该方法得到的市值规模因子体现出市值规模小的投资组合与市值规模大的投资组合之间的收益率差异，剔除账面市值比因素所造成的影响。

HMLt=BHt+SHt2-(BLt+SLt)2

同理，得到账面市值比因子也保证了解释变量只考虑账面市值比所产生的影响，反映账面市值比高的投资组合与账面市值比低的投资组合之间的收益率差异。

尽管许多西方研究证明FF3模型在美股已经失效，后续Fama本人也提出有效度更高的FF5模型，但是对于FF3的建模和选因子方式仍然具有借鉴意义。本文将尽可能呈现对FF3的理解的建模。

数据来源和数据清洗

本文数据来源为CSMAR数据库、网络完整的整理数据、中国资产管理研究中心的FF5的因子数据，包括每日收盘价、账面市值比、流通市值、无风险利率。时间区间为2017年1月1日至2020年12月31日四年时间的数据，频率为日，数据对象为中国A股所有股票，总计3391010行截面数据。

本文使用到的python包包括Pandas、Numpy、statsmodel、matplotlib、tqdm等包。

考虑到数据的有限来源以及无法获取tushare的pro接口，本文对于数据的清洗比较简陋，对于缺失值只是做一些简单填充和删除缺失值，并未考虑合适的插补数据方法。同时，考虑到股票的退市，理应对退市股票进行标记，并得到其退市区间，但由于数据过于简陋，且从CSMAR数据库导入该部分数据并不熟悉，最后选择放弃标记ST区间，仅仅对数据简单清洗，侧重于对FF3的建模和计算筛选股票alpha并可视化结果。

按照FF3的方法，本文需要对所有个股数据进行标记，分别按照每年五月末的市值分为两组，前50%是Big组，后50%是Small组；同时按照BM账面市值比分为三组，前30%是High组，后30%是Low组，其余的都是Medium组；一旦分组形成，以年为周期不发生改变，直到下一年末再重新分组。但是，我国上市公司企业年报的披露一般都在4月末，因此在构建因子时本文选用5月初的数据来分组。得到分组后的标记数据后，与所有股票数据进行合并。

随后，根据合并数据计算每个股票的收益率。同时，将每个月按照SIZE与BM分组，计算组内各支股票收益率市值加权综合收益率。

尽管由于经济形势的变化，普遍都认为FF3的因子建模得到的效果并不好，但是如果默认三因子模型的有效性，且市场风险、市值风险、账面市值比，这三类风险能很好地解释个股的超额收益，alphai的长期均值应该是0。那么，如果对于某个时期的股票，回归得到alphai<0，说明这段时间里面收益率偏低（因此股价也偏低），而根据有效市场假设，偏离值总会回归到均值。也就是说，如果在某个时刻根据FF3计算建模出的alpha值进行选股，选出alpha最小的N支股票，持仓一段时间后，下一个调仓日重复选股计算alpha最小的N支股票并持有，就可能有超额收益。

然而，由于我没有聚宽和其他量化的接口，也并未掌握有效的回测框架，导致学习曲线异常陡峭，于是放弃这种63日持仓择股的方法，打算后续攒积分后逐步开发这些策略。因此，我并没有按照选股是按如此进行的，为了方便循环建模，只是单纯计算如果按照alpha最大的股票进行计算收益率，持仓期设置为1年（也就是FF3计算因子分组的一年），每年一循环，找到前二十个alpha最大的股票，最终组合的收益率的结果会是多少。本着一种死马当活马医的想法，对所有时间节点的股票进行建模。

主要代码

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import tushare as ts
from pylab import mpl
import statsmodels.api as sm
import statsmodels.stats.api as sms
from tqdm import tqdm
import pickle
mpl.rcParams['font.sans-serif']=['SimHei']
mpl.rcParams['axes.unicode_minus']=False

#Fama三因子策略
# bm.csv 账面市值比
# cap.csv 流通市值
# close.csv 收盘价

#读取csv格式数据
bm = pd.read_csv("bm.csv")
cap = pd.read_csv("cap.csv")
close = pd.read_csv("close.csv")
#将data一列改为时间索引
bm = bm.set_index("Date")
cap = cap.set_index('Date')
close = close.set_index('Date')
#将所有数据均按照date格式输出
datelist=pd.to_datetime(cap.index.astype("int").astype("str"))
close = close.set_index(datelist)
bm = bm.set_index(datelist)
cap = cap.set_index(datelist)

#使用stack函数将数据变为标准面板数据并将三个dataframe并到一起
close_1 = close.stack().reset_index(name='close') #重新设置ind索引
bm_1 = bm.stack().reset_index(name='B/M')
cap_1 = cap.stack().reset_index(name='mkt')
close_1.columns = ['Date','stock_name','close']
bm_1.columns = ['Date','stock_name','B/M']
close_1 = close_1.join(bm_1, lsuffix='_close',rsuffix='_bm') #使用join函数合并列
close_1 = close_1.join(cap_1,lsuffix='_close',rsuffix='_cap')
my_data = close_1.drop(columns=['Date_bm','stock_name_bm','level_1','Date'])#删掉多余的列
my_data.columns = ['Date','code','close','B/M','mkt'] #得到初始数据

#去ST对数据要求很高，故直接填充
my_data = my_data.fillna(method='ffill', axis=0, inplace=False)
my_data = my_data.set_index('Date')
my_data

#定义两个画图的函数方便处理
def bar_plot_describe(df):
    df.plot(kind='bar', legend=True, figsize=(12, 8))
    plt.title('Results')
    plt.ylabel('Value')
    plt.show()

def box_plot_describe(df):
    df.boxplot(vert=False, grid=True, figsize=(12, 8))
    plt.title('Results')
    plt.show()

#创建一个新列包含日期的年份和月份
my_data['YearMonth'] = my_data.index.to_period('M')

#通过找到每个月第一个交易日的索引来确认
first_day_index = my_data[my_data['YearMonth'].dt.month == 5].groupby('YearMonth').apply(lambda x: x.index.min())
may_start_data = my_data.loc[first_day_index]

#得到包含五月初的所有数据
may_start_data =may_start_data.reset_index()

#分组
def mark_FF_group(df):
    df = df[df['close'] != 0]
    df['Date'] = pd.to_datetime(df['Date'])
    
    # 根据年份进行分组
    df['Year'] = df['Date'].dt.year

    # 根据市值排序，分为Big和Small组
    df['Mkt_Rank'] = df.groupby('Year')['mkt'].rank(pct=True)
    df['Size_Group'] = 'S'
    df.loc[df['Mkt_Rank'] >= 0.5, 'Size_Group'] = 'B'

    # 根据BM账面市值比排序，分为High、Medium、Low组
    df['BM_Rank'] = df.groupby('Year')['B/M'].rank(pct=True)
    df['BM_Group'] = 'M'
    df.loc[df['BM_Rank'] <= 0.3, 'BM_Group'] = 'L'
    df.loc[df['BM_Rank'] > 0.7, 'BM_Group'] = 'H'
    
    #同时标记
    df['portflio_mark'] =df['Size_Group'] +'/'+df['BM_Group']
    return df

marked_df = mark_FF_group(may_start_data)
marked_df['portflio_date'] = pd.to_datetime(marked_df['Date'],format="%Y-%m-%d").dt.strftime("%Y%m%d").astype(int).apply(lambda x: x//100)

#将标记好的数据与股票数据合并
my_data = my_data.reset_index()
data_rt0 = pd.DataFrame(my_data,columns=['Date','code','close','B/M','mkt'])
data_rt0['Date'] = pd.to_datetime(data_rt0['Date'],format="%Y-%m-%d").dt.strftime("%Y%m%d").astype(int)
data_rt0['portflio_date'] = data_rt0['Date'].apply(lambda x:x//10000*100+5 if (x//100%10>=5) or (x//1000%10==1) else (x//10000-1)*100+5)
data_rt1 = pd.merge(data_rt0,marked_df.drop(['Date','B/M','mkt','close'],axis=1),how='inner',left_on=['code','portflio_date'],right_on=['code','portflio_date'])
data_rt1

#计算得到每个股票的收益率序列
data_rt = data_rt1.sort_values(by='Date',ascending=True)
data_rt['ret'] = data_rt.groupby('code')['close'].apply(lambda x:x.pct_change()).dropna(axis =0)
data_rt.dropna(axis =0)

#将每个月按照SIZE与BM分组，计算组内各支股票收益率市值加权综合收益率ret
port_ret = data_rt.groupby(['Date','portflio_mark']).apply(lambda x: (x['ret']*x['mkt']).sum()/x['mkt'].sum())
port_ret = port_ret.reset_index()
port_ret.rename(columns={port_ret.columns[-1]:'ret'},inplace=True)  #重新设置列名
port_ret = port_ret.pivot(index='Date',columns='portflio_mark',values='ret')  #将表格改为透视表
port_ret

#计算三因子的SMB和HML
port_ret['SMB']=(port_ret['S/L']+port_ret['S/M']+port_ret['S/H'])/3-(port_ret['B/L']+port_ret['B/M']+port_ret['B/H'])/3
port_ret['HML']=(port_ret['S/H']+port_ret['B/H'])/2-(port_ret['S/L']+port_ret['B/L'])/2
ff3=port_ret.loc[:,['SMB','HML']].copy()
ff3=ff3.reset_index()
ff3.columns=['DATE','SMB','HML']
ff3['DATE'] = pd.to_datetime(ff3['DATE'].astype(str),format='%Y%m%d')
ff3     

#为了方便计算，fama三因子的市场风险溢价因子由其他导入
fama_5 = pd.read_csv("D:/桌面/python数据分析/Fama五因子数据/中国资产管理研究中心数据/fivefactor_daily.csv",index_col='trddy',parse_dates=True)
fama_5 = fama_5.loc['2017-05-02':'2020-12-31']
fama_5['portflio_date'] = pd.to_datetime(fama_5.index,format='%Y-%m-%d').strftime("%Y%m%d").astype(int)
fama_5 = fama_5.reset_index()
fama_5 = fama_5.rename(columns={'trddy': 'DATE'})
ff3 = pd.merge(ff3, fama_5[['DATE', 'mkt_rf','rf']], on='DATE', how='left')

# 绘制三因子的收益率折线图
df_resampled = ff3[['SMB', 'HML', 'mkt_rf', 'rf']].cumsum()
plt.figure(figsize=(12, 10))
plt.plot(df_resampled.index, df_resampled['SMB'], label='SMB', marker='o', linestyle='-',markevery=5)
plt.plot(df_resampled.index, df_resampled['HML'], label='HML', color ='lightblue',marker='o', linestyle='-',markevery=5)
plt.plot(df_resampled.index, df_resampled['mkt_rf'], label='MKT', color= 'red',marker='o', linestyle='-',markevery=5)
plt.plot(df_resampled.index, df_resampled['rf'], label='Rf', color='green',marker='o', linestyle='-',markevery=5)
plt.xlabel('DATE')
plt.ylabel('Values')
plt.title('Time Series Plot')
plt.legend()
plt.tight_layout()
plt.show()

new_df = data_rt[['Date','code','ret',]]
new_df = new_df.rename(columns = {'Date':'DATE'})
new_df['DATE'] = pd.to_datetime(new_df['DATE'].astype(int).astype(str),format='%Y%m%d')

# 得到合并后的所有数据
merged_df = pd.merge(new_df,ff3,on='DATE')
merged_df = merged_df.dropna(axis =0)
merged_df.isnull().sum() #判断是否有空值

# 创建几个空字典来保存回归结果
results_dict = {}
alpha_list = []
portfolio_returns = []
top_stocks = []
model_all = {}
unique_codes = merged_df['code'].unique()
merged_df.loc[:,'year'] = pd.DatetimeIndex(merged_df['DATE']).year
merged_df

#定义一个保存回归模型到本地的函数
def save_data(data,filename):
    with open(f'{filename}.pkl', 'wb') as file:
        pickle.dump(data, file)

#定义一个加载到本地的这些文件的函数
def load_data(filename):
    with open(f'{filename}.pkl', 'rb') as file:
        data = pickle.load(file)
    return data

#每年一循环，找到前二十个alpha最大的股票
for year in tqdm(merged_df['year'].unique(), desc='Processing years'):
    alpha_list = []
    for code in tqdm(unique_codes,desc ='Processing stocks'):
        try: #try很有用！如果成功就执行，不成功就except兜底！
            data = merged_df[(merged_df['code'] == code) & (merged_df['year'] == year)]
            data = data.dropna(axis =0)

            X = data[['SMB', 'HML', 'mkt_rf']]
            X = sm.add_constant(X)  # 添加截距项
            y = data['ret']
            #担心数据量，因此提前筛选一些股票
            if len(y) < 4 or len(X) < 4:
                continue

            # 回归！
            model = sm.OLS(y, X, missing='drop')
            results = model.fit()
            alpha = results.params['const']
            alpha_list.append({'code': code, 'alpha': alpha})
            
            # 将回归结果保存到字典中
            results_dict[(code,year)] = results.params
            model_all[(code,year)] = results
        except Exception as e:
            print(f"Error occurred while processing stock {code} in year {year}: {e}")
            save_data(results_dict, 'results_dict_backup')
            continue
    #得到一些数据
    top_stock20 = sorted(alpha_list, key=lambda x: x['alpha'], reverse=True)[:20]
    top_stocks.append({'year':year,'top_stock20':top_stock20})
    portfolio_return = np.mean([stock['alpha'] for stock in top_stock20])
    portfolio_returns.append({'year': year, 'return': portfolio_return})

save_data(results_dict, 'results_dict')
save_data(portfolio_returns, 'portfolio_returns')
save_data(top_stocks, 'top_stocks')
save_data(model_all,'model_all')

#提取这些数据，并接下来分开查看方便可视化
results_dict = load_data('results_dict')
portfolio_returns = load_data('portfolio_returns')
top_stocks = load_data('top_stocks')
model_all = load_data('model_all')

rsquared_list = [{'code': code, 'year': year, 'rsquared': result.rsquared} for (code, year), result in model_all.items()]
top_5_stocks = sorted(rsquared_list, key=lambda x: x['rsquared'], reverse=True)[:10]
for stock in top_10_stocks:
    print(f"Stock: {stock['code']}, Year: {stock['year']}, R-squared: {stock['rsquared']}")