【问题标题】:Fama Macbeth Regression in Python (Pandas or Statsmodels)Python 中的 Fama Macbeth 回归(Pandas 或 Statsmodels)
【发布时间】:2014-07-27 06:57:41
【问题描述】:

计量经济学背景

Fama Macbeth 回归是指对面板数据运行回归的过程(其中有 N 个不同的个体,每个个体对应于多个时期 T,例如日、月、年)。所以总共有 N x T obs。请注意,如果面板数据不平衡,则可以。

Fama Macbeth 回归是首先跨部门运行每个时期的回归,即将给定时期 t 内的 N 个个体汇集在一起​​。并为 t=1,...T 执行此操作。所以总共运行了 T 个回归。然后我们有每个自变量的时间序列系数。然后我们可以使用系数的时间序列进行假设检验。通常我们取平均值作为每个自变量的最终系数。我们使用 t-stats 来检验显着性。

我的问题

我的问题是在熊猫中实现这一点。从 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索引的数据框,并且数据框的每一列都应该包含每个变量interceptvar1var2var3的系数。

我也检查了statsmodels,他们也没有这样的内置程序。

是否有任何软件包可以生成出版质量的回归表?像 Stata 中的 outreg2 和 R 中的 texreg? 谢谢你的帮助!

【问题讨论】:

  • 到问题的最后一部分:statsmodels 有一个summary_col 函数,但它缺少一些选项(和花里胡哨)github.com/statsmodels/statsmodels/issues/1637
  • 谢谢!但它看起来很粗糙,而且不像 R 的 texreg 或 stata 的 outreg2 那样整洁。我想我将不得不切换到 R 输出

标签: python r pandas statsmodels


【解决方案1】:

更新以反映 Fama-MacBeth 截至 2018 年秋季的图书馆情况。fama_macbeth 函数已从 pandas 中删除一段时间。那么你有什么选择呢?

  1. 如果您使用的是 python 3,那么您可以在 LinearModels 中使用 Fama-MacBeth 方法:https://github.com/bashtage/linearmodels/blob/master/linearmodels/panel/model.py

  2. 如果您使用的是 python 2 或者只是不想使用 LinearModels,那么最好的选择可能是自己动手。

例如,假设您在如下面板中有 Fama-French 行业投资组合(您还计算了一些变量,例如过去的 beta 或过去的回报率以用作您的 x 变量):

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结果对象方法@987654330 @ 没有返回 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。编写一个新函数来进行 ols 估计...类似于以下内容(请注意,我没有做任何检查这些矩阵的等级之类的事情...):

def ols_np(data,yvar,xvar):
    gamma,_,_,_ = np.linalg.lstsq(data[xvar],data[yvar],rcond=None)
    return pd.Series(gamma)

如果您仍在使用旧版本的 pandas,以下方法将起作用:

下面是在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-var 和 x-vars 具有多索引,其中日期作为第一个变量,股票/公司/实体 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-var 必须是 dataframe,因此如果您只传递一列,则需要将其传递为 df[['x']]

如果您不想拦截,您必须这样做:

>>> fm  = pd.fama_macbeth(y=df['y'],x=df[['x']],intercept=False)

【讨论】:

  • 对不起,我的意思是把它移到statsmodels
  • 您是否愿意在groupby 中告诉我如何手动执行此操作,即我在自己的代码的第二部分中犯了什么错误?
  • 我的统计模型是最新的。我想将每个系数保持为一行中的不同列,应该由date 索引。我会尝试添加pd.Series 看看是否有效
  • 您的代码适用于我的简单示例:res = df.groupby('date').apply(fmreg,'y~1+x')。当我添加另一个 x 变量时它也可以正常工作:df.groupby('date').apply(fmreg,'y~x + x2')
  • 对,我得到了你想要的输出,而没有真正修改你的代码。除非您展示一个可重现的小问题示例,否则我不确定我是否真的能提供帮助。
【解决方案2】:

编辑:新图书馆

存在一个更新的库,可以通过以下命令安装:

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)

注意:关闭 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

这些表格可以导出到 LaTeX 以便于写入:

tbl.to_latex()

【讨论】:

  • 这些byu链接还存在吗?
  • @Martien 我刚刚更新了链接和安装信息。文档移至 readthedocs 上的新位置,并从测试 PyPi 迁移到普通 PyPI。
【解决方案3】:

一种快速而肮脏的解决方案,可以解决问题并继续使用您正在使用的相同东西。

它对我有用。

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

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

【讨论】:

    猜你喜欢
    • 2020-08-10
    • 2012-04-19
    • 1970-01-01
    • 1970-01-01
    • 2015-07-22
    • 2016-01-24
    • 2015-11-05
    • 2017-08-17
    • 2016-09-07
    相关资源
    最近更新 更多