如果你经常用stata写论文,会了解stata有个outreg2的函数,可以把回归的结果输出成非常规范的论文格式,并且可以把多个回归结果并在一起,方便对比。例如下图
本文的目的是用python实现outreg2的效果,得到上面这样的结果,方便对比和分析。
01
源码说明
其实也不用自己手动写,statsmodels模块里有一个summary_col函数,可以实现以上的功能,不过效果没有stata那么好,毕竟python也不是专业的计量分析软件,但好在代码并不难,所以如果有一些个性化的需求,自己改一改也挺容易的。
先解释代码,再上例子。首先看看summary_col的说明:
一共有七个参数,一一解释(源码也不难,有兴趣可以自己看看)。
01
OLS格式化输出
例子使用前文Fama-Macbeth中使用过的数据,首先取其中一期的数据做回归,这里主要是展示格式化输出的结果,所以不要太在意系数的符号和显著性。
读入数据
price = pd.read_csv('price.csv',index_col = 0)
pb = pd.read_csv('pb.csv',index_col = 0)
roe = pd.read_csv('roe.csv',index_col = 0)
mkt = pd.read_csv('mkt.csv',index_col = 0)
ST = pd.read_excel('ST.xlsx')
ind = pd.read_csv('中信一级行业.csv',encoding = 'gbk')
'''
收益率计算
'''
price['tradedate'] = price.tradedate.apply(getdate)
pb['tradedate'] = pb.tradedate.apply(getdate)
roe['tradedate'] = roe.tradedate.apply(getdate)
mkt['tradedate'] = mkt.tradedate.apply(getdate)
ind['tradedate'] = ind.tradedate.apply(getdate)
ret_m = getRet(price,freq ='m',if_shift = True)
mom1 = getRet(price,freq ='m',if_shift = False)
mom1 = mom1.rename(columns = {'ret':'mom1'})
fall = pd.merge(pb,roe,left_on = ['tradedate','stockcode'],right_on = ['tradedate','stockcode'])
fall = pd.merge(fall,mkt,left_on = ['tradedate','stockcode'],right_on = ['tradedate','stockcode'])
fall['tradedate'] = fall.tradedate.apply(getdate)
fall = pd.merge(fall,mom1,left_on = ['tradedate','stockcode'],right_on = ['tradedate','stockcode'])
del fall['rptdate']
ind,on = ['tradedate','stockcode'])
alldata = pd.merge(fall,ret_m,left_on = ['tradedate','stockcode'],right_on = ['tradedate','stockcode'])
alldata['mktcap'] = np.log(alldata.mktcap)
# classname转虚拟变量
alldata = pd.get_dummies(alldata,columns = ['classname'],drop_first = True,prefix = '',prefix_sep = '')
indname = list(alldata.columns[-28:])
from statsmodels.iolib.summary2 import summary_col
取一个截面上的数据,分别做五次回归:
最后把五次回归的结果合并在一起格式化输出,注意这里行业用的时中信一级行业,虚拟变量个数很多,所以用drop_omitted设置不输出这些虚拟变量的系数。
代码如下
#### OLS ####
fmdata = alldata.set_index(['stockcode','tradedate'])
fmdata = fmdata.fillna(0)
ols_data = fmdata.xs(datetime.date(2010,1,29),level = 1)
from statsmodels.iolib.summary2 import *
y = ols_data['ret']
x1 = ols_data[['pb'] + indname]
x2 = ols_data[['mktcap'] + indname]
x3 = ols_data[['mom1'] + indname]
x4 = ols_data[['roe_ttm'] + indname]
x5 = ols_data[['pb','mktcap','mom1','roe_ttm'] + indname]
res_ols1 = sm.OLS(y,sm.add_constant(x1)).fit()
res_ols2 = sm.OLS(y,sm.add_constant(x2)).fit()
res_ols3 = sm.OLS(y,sm.add_constant(x3)).fit()
res_ols4 = sm.OLS(y,sm.add_constant(x4)).fit()
res_ols5 = sm.OLS(y,sm.add_constant(x5)).fit()
res_ols1.summary()
summary_col([res_ols1,res_ols2,res_ols3,res_ols4,res_ols5],
model_names = ['f1','f2','f3','f4','f5'],
stars = True,regressor_order = ['const','pb','mktcap','mom1','roe_ttm'],
drop_omitted = indname,info_dict = {'':lambda x: '',
'':lambda x: '',
'Observation':lambda x:str(int(x.nobs)),
})
结果如下
这里的info_dict里子定义了三个行数,效果是两行空白,第三行输出变量个数,也就是图里的Observation,如果你想在结果里输出更多的统计量,也可以用类似的方法实现,小数位数也是可以调的。
02
Fama-Macbeth 格式化输出
Fama-macbeth之前提到可以使用FamaMacBeth函数实现,但如果直接用上面的方式对Fama-macbeth的结果输出会有一些问题。
首先还是做上面五次回归
y = fmdata['ret']
x1 = fmdata[['pb'] + indname]
x2 = fmdata[['mktcap'] + indname]
x3 = fmdata[['mom1'] + indname]
x4 = fmdata[['roe_ttm'] + indname]
x5 = fmdata[['pb','mktcap','mom1','roe_ttm'] + indname]
res_fm1 = FamaMacBeth(y,sm.add_constant(x1)).fit(cov_type= 'kernel',debiased = False,bandwidth = 5)
res_fm2 = FamaMacBeth(y,sm.add_constant(x2)).fit(cov_type= 'kernel',debiased = False,bandwidth = 5)
res_fm3 = FamaMacBeth(y,sm.add_constant(x3)).fit(cov_type= 'kernel',debiased = False,bandwidth = 5)
res_fm4 = FamaMacBeth(y,sm.add_constant(x4)).fit(cov_type= 'kernel',debiased = False,bandwidth = 5)
res_fm5 = FamaMacBeth(y,sm.add_constant(x5)).fit(cov_type= 'kernel',debiased = False,bandwidth = 5)
如果这时候使用和上面同样的方式去输出会发现报错了
这个去看看源码会发现是ols的属性里有bse,fama-macbeth的属性里没有bse,但fm也有同样的统计量,只是名称不一样,所以这里需要多加一步转换函数来对fm的回归结果做一些转换,然后就可以实现输出了。
def transfor(resfm):
resfm.bse = resfm.std_errors
resfm.tvalues = resfm.tstats
resfm.model.exog_names = list(resfm.model.exog.dataframe.columns)
resfm.model.endog_names = list(resfm.model.dependent.dataframe.columns)
return resfm
res_fm1 = transfor(res_fm1)
res_fm2 = transfor(res_fm2)
res_fm3 = transfor(res_fm3)
res_fm4 = transfor(res_fm4)
res_fm5 = transfor(res_fm5)
summary_col([res_fm1,res_fm2,res_fm3,res_fm4,res_fm5],
model_names = ['f1','f2','f3','f4','f5'],
regressor_order = ['const','pb','mktcap','mom1','roe_ttm'],
stars = True,
drop_omitted = indname,info_dict = {'':lambda x: '',
'':lambda x: '',
'Observation':lambda x:str(int(x.nobs)),
})
最终输出结果如下
这里会发现R-squared Adj输出是nan,主要是fama-macbeth回归没有调整R2方的概念,可以自己设置不输出R2或者换成别的统计量。
另外这个包目前还是在完善过程中,所以如果python版本不一样,输出结果可能会有一些差异,比如上图是用python3.7实现的,python3.8实现出来R2的结果会显示在回归系数的下方。当然这些都是小问题了,可以根据自己的需要去调整。