Python中的Fama Macbeth回归（Pandas或Statsmodels）

Question

经济计量背景

Fama-Macbeth回归是指对面板数据进行回归的过程（其中有N个不同的个体，每个个体对应多个时段T，例如天、月、年）。所以总的来说有N x T obs。请注意，如果面板数据不平衡，则可以。

Fama-Macbeth回归是对每个时期进行跨部门的第一次回归，即在给定的时期t内将N个个体聚集在一起，并在t=1，…t时进行回归，因此在总t回归中进行回归。然后我们得到每个自变量的系数时间序列。然后利用系数的时间序列进行假设检验。通常我们取平均值作为每个自变量的最终系数。我们用t统计量来检验显著性。

我的问题

我的问题是在熊猫身上实施。从pandas的源代码中，我注意到有一个名为fama_macbeth的过程。但是我找不到任何关于这个的文件。

这个操作也可以通过groupby轻松完成。目前我正在做：

def fmreg(data,formula):
    return smf.ols(formula,data=data).fit().params[1]

res=df.groupby('date').apply(fmreg,'ret~var1')

这是这样的，res是一个由date索引的序列，序列的值是params[1]，这是var1的系数。但现在我想有更多的自变量，我需要提取所有这些自变量的系数，但我无法计算出来。我试过这个

def fmreg(data,formula):
    return smf.ols(formula,data=data).fit().params

res=df.groupby('date').apply(fmreg,'ret~var1+var2+var3')

这行不通。期望的结果是res是由date索引的数据帧，数据帧的每一列应该包含每个变量intercept、var1、var2和var3的系数。

我还用statsmodels检查过，它们也没有这样的内置过程。

是否有任何包可以生成发布质量回归表？就像Stata中的outreg2和R中的texreg？ 谢谢你的帮助！

Answer 1

编辑：新建库

已存在可通过以下命令安装的更新库：

pip install finance-byu

这里的文档：https://fin-585-byu.readthedocs.io/en/latest/contents.html

新的库包括Fama Macbeth回归实现，速度得到了提高，并且更新了Regtable类。新的图书馆还包括GRS统计，这也可能是有用的人对Fama麦克白回归感兴趣。

---

作为对Karl D.上述答案的更新，现在有一个非常年轻的库，它在python中生成回归表，并包含Fama Macbeth回归函数。

文档：https://byu-finance-library-finance-byu.readthedocs.io/en/latest/contents.html

当前可以从测试PyPi服务器导入库：

pip install -i https://test.pypi.org/simple/ finance-byu

文档中的这一页概述了Fama Macbeth函数：https://byu-finance-library-finance-byu.readthedocs.io/en/latest/fama_macbeth.html

有一个实现非常类似于Karl D.在上面使用numpy的线性代数函数的实现，该实现利用joblib进行并行化，以在数据中有大量时间段时提高性能，以及一个使用numba进行优化的实现，它在小数据集上削减了一个数量级。

下面是一个小的模拟数据集示例，如文档中所示：

>>> from finance_byu.fama_macbeth import fama_macbeth, fama_macbeth_parallel, fm_summary, fama_macbeth_numba
>>> import pandas as pd
>>> import time
>>> import numpy as np
>>> 
>>> n_jobs = 5
>>> n_firms = 1.0e2
>>> n_periods = 1.0e2
>>> 
>>> def firm(fid):
>>>     f = np.random.random((int(n_periods),4))
>>>     f = pd.DataFrame(f)
>>>     f['period'] = f.index
>>>     f['firmid'] = fid
>>>     return f
>>> df = [firm(i) for i in range(int(n_firms))]
>>> df = pd.concat(df).rename(columns={0:'ret',1:'exmkt',2:'smb',3:'hml'})
>>> df.head()

        ret     exmkt       smb       hml  period  firmid
0  0.766593  0.002390  0.496230  0.992345       0       0
1  0.346250  0.509880  0.083644  0.732374       1       0
2  0.787731  0.204211  0.705075  0.313182       2       0
3  0.904969  0.338722  0.437298  0.669285       3       0
4  0.121908  0.827623  0.319610  0.455530       4       0

>>> result = fama_macbeth(df,'period','ret',['exmkt','smb','hml'],intercept=True)
>>> result.head()

        intercept     exmkt       smb       hml
period                                         
0        0.655784 -0.160938 -0.109336  0.028015
1        0.455177  0.033941  0.085344  0.013814
2        0.410705 -0.084130  0.218568  0.016897
3        0.410537  0.010719  0.208912  0.001029
4        0.439061  0.046104 -0.084381  0.199775

>>> fm_summary(result)

               mean  std_error      tstat
intercept  0.506834   0.008793  57.643021
exmkt      0.004750   0.009828   0.483269
smb       -0.012702   0.010842  -1.171530
hml        0.004276   0.010530   0.406119

>>> %timeit fama_macbeth(df,'period','ret',['exmkt','smb','hml'],intercept=True)
123 ms ± 117 µs per loop (mean ± std. dev. of 7 runs, 10 loops each  
>>> %timeit fama_macbeth_parallel(df,'period','ret',['exmkt','smb','hml'],intercept=True,n_jobs=n_jobs,memmap=False)  
146 ms ± 16.9 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
>>> %timeit fama_macbeth_numba(df,'period','ret',['exmkt','smb','hml'],intercept=True)
5.04 ms ± 5.2 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

注意：关闭memmap可以进行公平比较，而不必每次运行时生成新数据。使用memmap，并行实现只需提取缓存的结果。

下面是表类的两个简单实现，它们也使用模拟数据：

>>> from finance_byu.regtables import Regtable
>>> import pandas as pd
>>> import statsmodels.formula.api as smf
>>> import numpy as np
>>> 
>>> 
>>> nobs = 1000
>>> df = pd.DataFrame(np.random.random((nobs,3))).rename(columns={0:'age',1:'bmi',2:'hincome'})
>>> df['age'] = df['age']*100
>>> df['bmi'] = df['bmi']*30
>>> df['hincome'] = df['hincome']*100000
>>> df['hincome'] = pd.qcut(df['hincome'],16,labels=False)
>>> df['rich'] = df['hincome'] > 13
>>> df['gender'] = np.random.choice(['M','F'],nobs)
>>> df['race'] = np.random.choice(['W','B','H','O'],nobs)
>>> 
>>> regformulas =  ['bmi ~ age',
>>>                 'bmi ~ np.log(age)',
>>>                 'bmi ~ C(gender) + np.log(age)',
>>>                 'bmi ~ C(gender) + C(race) + np.log(age)',
>>>                 'bmi ~ C(gender) + rich + C(gender)*rich + C(race) + np.log(age)',
>>>                 'bmi ~ -1 + np.log(age)',
>>>                 'bmi ~ -1 + C(race) + np.log(age)']
>>> reg = [smf.ols(f,df).fit() for f in regformulas]
>>> tbl = Regtable(reg)
>>> tbl.render()

产生以下结果：

>>> df2 = pd.DataFrame(np.random.random((nobs,10)))
>>> df2.columns = ['t0_vw','t4_vw','et_vw','t0_ew','t4_ew','et_ew','mktrf','smb','hml','umd']
>>> regformulas2 = ['t0_vw ~ mktrf',
>>>                't0_vw ~ mktrf + smb + hml',
>>>                't0_vw ~ mktrf + smb + hml + umd',
>>>                't4_vw ~ mktrf',
>>>                't4_vw ~ mktrf + smb + hml',
>>>                't4_vw ~ mktrf + smb + hml + umd',
>>>                'et_vw ~ mktrf',
>>>                'et_vw ~ mktrf + smb + hml',
>>>                'et_vw ~ mktrf + smb + hml + umd',
>>>                't0_ew ~ mktrf',
>>>                't0_ew ~ mktrf + smb + hml',
>>>                't0_ew ~ mktrf + smb + hml + umd',
>>>                't4_ew ~ mktrf',
>>>                't4_ew ~ mktrf + smb + hml',
>>>                't4_ew ~ mktrf + smb + hml + umd',
>>>                'et_ew ~ mktrf',
>>>                'et_ew ~ mktrf + smb + hml',
>>>                'et_ew ~ mktrf + smb + hml + umd'
>>>                ]
>>> regnames = ['Small VW','','',
>>>             'Large VW','','',
>>>             'Spread VW','','',
>>>             'Small EW','','',
>>>             'Large EW','','',
>>>             'Spread EW','',''
>>>             ]
>>> reg2 = [smf.ols(f,df2).fit() for f in regformulas2]
>>> 
>>> tbl2 = Regtable(reg2,orientation='horizontal',regnames=regnames,sig='coeff',intercept_name='alpha',nobs=False,rsq=False,stat='se')
>>> tbl2.render()

产生以下结果：

Regtable类的文档在这里：https://byu-finance-library-finance-byu.readthedocs.io/en/latest/regtables.html

这些表格可以导出为乳胶，便于写入：

tbl.to_latex()