[解決済み] PythonによるFama Macbeth回帰 (PandasまたはStatsmodels)

質問

経済学的背景

Fama Macbeth回帰とは、パネルデータ（N人の異なる個人がいて、各個人が複数の期間T、例えば、日、月、年に対応する）に対して回帰を実行する手順を指します。つまり、合計でN x T個のオブザーバが存在します。パネルデータがバランスしていなくても問題ないことに注意してください。

Fama Macbeth回帰は、まず各期間について横断的に回帰を実行する、つまり、ある期間tにN人の個人をプールしておく。そして、これをt=1,...Tについて行う。したがって、合計T回の回帰が実行される。そして、各独立変数の係数の時系列が得られます。そして、係数の時系列を用いて仮説検定を行うことができます。通常、我々は、各独立変数の最終的な係数として、平均をとります。そして、t-statsで有意性を検定します。

私の問題

私の問題は、これを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の値は 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 そのようなビルトインプロシージャはありません。

また、出版品質の回帰表を作成できるパッケージはないのでしょうか？例えば outreg2 はStataで、そして texreg をRで使うことはできますか？ よろしくお願いします。

どのように解決するのですか？

2018年秋時点のFama-MacBethのライブラリ状況を反映したアップデートです。その 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の産業別ポートフォリオがあるとします（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 結果オブジェクトメソッド 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 . olsの推定を行う新しい関数を書きます ... 以下のようなものです（これらの行列のランクをチェックするようなことは何もしていないことに注意してください ...):

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

また、まだ古いバージョンの pandas は、次のようにすれば動作します。

を使用する例です。 fama_macbeth 関数は pandas :

>>> 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-var が、最初の変数として日付、インデックスの 2 番目の変数として株式/企業/団体 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 として渡さなければならないので、1列だけ渡す場合は df[['x']] .

インターセプトがいらないなら、しなければならないこと。

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