用python输出stata一样的标准化回归结果

量化小白

发布于 2020-11-03 15:31:32

4.7K0

发布于 2020-11-03 15:31:32

文章被收录于专栏：量化小白上分记量化小白上分记

如果你经常用stata写论文，会了解stata有个outreg2的函数，可以把回归的结果输出成非常规范的论文格式，并且可以把多个回归结果并在一起，方便对比。例如下图

本文的目的是用python实现outreg2的效果，得到上面这样的结果，方便对比和分析。

源码说明

其实也不用自己手动写，statsmodels模块里有一个summary_col函数，可以实现以上的功能，不过效果没有stata那么好，毕竟python也不是专业的计量分析软件，但好在代码并不难，所以如果有一些个性化的需求，自己改一改也挺容易的。

先解释代码，再上例子。首先看看summary_col的说明：

一共有七个参数，一一解释（源码也不难，有兴趣可以自己看看）。

results：回归模型的结果，statsmodel回归的结果，多个结果用list装在一起；
float_format：数据显示的格式，默认四位小数；
stars：是否有*，True为有；
model_names：模型的名称，用字符串装在list里就可以了，输出时的表头，也可以不设置；
info_dict：这个用来生成一些自己想加在格式化输出里的内容，dict格式；
regressor_order：用来设置自变量的显示顺序；
drop_omitted：是否包含没有包括在regressor_order中的变量。

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

取一个截面上的数据，分别做五次回归：

ret~pb+行业虚拟变量
ret~mkt+行业虚拟变量
ret~mom1+行业虚拟变量
ret~roe_ttm+行业虚拟变量
ret~pb+mkt+mom+roe_ttm+行业虚拟变量

最后把五次回归的结果合并在一起格式化输出，注意这里行业用的时中信一级行业，虚拟变量个数很多，所以用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，如果你想在结果里输出更多的统计量，也可以用类似的方法实现，小数位数也是可以调的。

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)),
                                                 })

最终输出结果如下