Python中的Fama Macbeth回归（Pandas或Statsmodels）
计量经济学背景
Fama Macbeth 回归是一种用于面板数据的回归方法。面板数据的意思是有 N 个不同的个体，每个个体对应多个时间段 T，比如天、月、年。总的来说，就是有 N x T 个观察值。即使面板数据不平衡也是可以的。
Fama Macbeth 回归的步骤是，首先在每个时间段内进行横截面回归，也就是说在某个特定的时间 t，把 N 个个体放在一起进行回归。然后对每个时间段 t=1,...,T 都这样做。这样总共会进行 T 次回归。接着，我们会得到每个自变量的系数时间序列。然后可以利用这些系数进行假设检验。通常，我们会取平均值作为每个自变量的最终系数，并用 t 统计量来检验显著性。
我的问题
我现在的问题是如何在 pandas 中实现这个方法。在 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 为索引的 Series，Series 的值是 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？谢谢你的帮助！
3 个回答
这是一个简单粗暴的解决办法，可以让你继续使用之前的东西。
对我来说，这个方法有效。
def fmreg(data,formula):
return smf.ols(formula,data=data).fit().params[:]
res = df.groupby('date').apply(fmreg,'ret~var1+var2+var3')
编辑：新库
现在有一个更新过的库，可以通过以下命令安装：
pip install finance-byu
文档在这里：https://fin-library.readthedocs.io/en/latest/
这个新库包含了Fama Macbeth回归的实现，还有一个叫Regtable的类，可以帮助你报告结果。
文档中的这一页介绍了Fama Macbeth的功能：https://fin-library.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)
注意：关闭内存映射可以让比较更公平，而不需要在每次运行时生成新数据。使用内存映射时，并行实现会直接调用缓存的结果。
这里有几个简单的表格类实现，也使用了模拟数据：
>>> 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
这些表格可以导出为LaTeX格式，方便你在写作中使用：
tbl.to_latex()
这是关于Fama-MacBeth库在2018年秋季的更新情况。fama_macbeth这个功能已经从pandas中移除一段时间了。那么你还有哪些选择呢？
如果你使用的是Python 3，那么可以在LinearModels中使用Fama-MacBeth方法：https://github.com/bashtage/linearmodels/blob/master/linearmodels/panel/model.py
如果你使用的是Python 2，或者不想使用LinearModels，那么你可能需要自己动手实现一个。
例如，假设你有Fama-French行业投资组合的数据，像下面这样（你还计算了一些变量，比如过去的beta值或过去的收益，用作你的自变量）：
In [1]: import pandas as pd
import numpy as np
import statsmodels.formula.api as smf
In [4]: df = pd.read_csv('industry.csv',parse_dates=['caldt'])
df.query("caldt == '1995-07-01'")
In [5]: Out[5]:
industry caldt ret beta r12to2 r36to13
18432 Aero 1995-07-01 6.26 0.9696 0.2755 0.3466
18433 Agric 1995-07-01 3.37 1.0412 0.1260 0.0581
18434 Autos 1995-07-01 2.42 1.0274 0.0293 0.2902
18435 Banks 1995-07-01 4.82 1.4985 0.1659 0.2951
Fama-MacBeth主要是每个月计算相同的横截面回归模型，所以你可以用groupby来实现。你可以创建一个函数，接收一个dataframe（这个数据来自groupby）和一个patsy公式；然后它会拟合模型并返回参数估计值。这里是一个简单的实现示例（注意这是原提问者几年前尝试做的事情……不太清楚当时为什么没成功，可能是因为当时statsmodels的结果对象方法params没有返回pandas的Series，所以返回值需要显式转换为Series……在当前版本的pandas 0.23.4中是可以正常工作的）：
def ols_coef(x,formula):
return smf.ols(formula,data=x).fit().params
In [9]: gamma = (df.groupby('caldt')
.apply(ols_coef,'ret ~ 1 + beta + r12to2 + r36to13'))
gamma.head()
In [10]: Out[10]:
Intercept beta r12to2 r36to13
caldt
1963-07-01 -1.497012 -0.765721 4.379128 -1.918083
1963-08-01 11.144169 -6.506291 5.961584 -2.598048
1963-09-01 -2.330966 -0.741550 10.508617 -4.377293
1963-10-01 0.441941 1.127567 5.478114 -2.057173
1963-11-01 3.380485 -4.792643 3.660940 -1.210426
然后计算均值、均值的标准误差和t检验（或者你想要的其他统计量）。像下面这样：
def fm_summary(p):
s = p.describe().T
s['std_error'] = s['std']/np.sqrt(s['count'])
s['tstat'] = s['mean']/s['std_error']
return s[['mean','std_error','tstat']]
In [12]: fm_summary(gamma)
Out[12]:
mean std_error tstat
Intercept 0.754904 0.177291 4.258000
beta -0.012176 0.202629 -0.060092
r12to2 1.794548 0.356069 5.039896
r36to13 0.237873 0.186680 1.274230
提高速度
使用statsmodels进行回归会有较大的开销（特别是你只需要估计的系数）。如果你想提高效率，可以从statsmodels切换到numpy.linalg.lstsq。写一个新的函数来进行最小二乘估计……像下面这样（注意我没有检查这些矩阵的秩……）：
def ols_np(data,yvar,xvar):
gamma,_,_,_ = np.linalg.lstsq(data[xvar],data[yvar],rcond=None)
return pd.Series(gamma)
如果你仍在使用旧版本的pandas，下面的代码将会有效：
这里是一个使用fama_macbeth函数的示例：
>>> df
y x
date id
2012-01-01 1 0.1 0.4
2 0.3 0.6
3 0.4 0.2
4 0.0 1.2
2012-02-01 1 0.2 0.7
2 0.4 0.5
3 0.2 0.1
4 0.1 0.0
2012-03-01 1 0.4 0.8
2 0.6 0.1
3 0.7 0.6
4 0.4 -0.1
注意结构。fama_macbeth函数期望y变量和x变量有一个多重索引，日期作为第一个变量，股票/公司/实体的ID作为索引的第二个变量：
>>> fm = pd.fama_macbeth(y=df['y'],x=df[['x']])
>>> fm
----------------------Summary of Fama-MacBeth Analysis-------------------------
Formula: Y ~ x + intercept
# betas : 3
----------------------Summary of Estimated Coefficients------------------------
Variable Beta Std Err t-stat CI 2.5% CI 97.5%
(x) -0.0227 0.1276 -0.18 -0.2728 0.2273
(intercept) 0.3531 0.0842 4.19 0.1881 0.5181
--------------------------------End of Summary---------------------------------
注意，直接打印fm会调用fm.summary
>>> fm.summary
----------------------Summary of Fama-MacBeth Analysis-------------------------
Formula: Y ~ x + intercept
# betas : 3
----------------------Summary of Estimated Coefficients------------------------
Variable Beta Std Err t-stat CI 2.5% CI 97.5%
(x) -0.0227 0.1276 -0.18 -0.2728 0.2273
(intercept) 0.3531 0.0842 4.19 0.1881 0.5181
--------------------------------End of Summary---------------------------------
另外，fama_macbeth函数会自动添加一个截距（与statsmodels的例程不同）。而且x变量必须是一个dataframe，所以如果你只传递一列，需要以df[['x']]的形式传递。
如果你不想要截距，你需要这样做：
>>> fm = pd.fama_macbeth(y=df['y'],x=df[['x']],intercept=False)