Python での Fama-Macbeth 回帰
- ファマ・マクベス回帰とその重要性
- Python で Fama-Macbeth 回帰を実装する手順
-
LinearModelsを使用して Python で Fama-Macbeth 回帰を実装する - Python で Fama-Macbeth 回帰を実装するための代替アプローチ
今日の記事では、Fama-Macbeth 回帰、その重要性、およびその実装について説明します。
ファマ・マクベス回帰とその重要性
資産価格理論では、リスク要因を使用して資産のリターンを説明します。 これらのリスク要因は、ミクロ経済的またはマクロ経済的である可能性があります。
ミクロ経済的なリスク要因には、企業規模と企業のさまざまな財務指標が含まれますが、マクロ経済的なリスク要因には、消費者のインフレと失業があります。
Fama-Macbeth 回帰は、資産価格モデルのテストに使用される 2 段階の回帰モデルです。 これらのリスク要因がポートフォリオや資産のリターンをどの程度正確に表しているかを測定するのは、実用的なアプローチです。
このモデルは、これらの要因へのエクスポージャーに関連するリスク プレミアムを決定するのに役立ちます。
要点は、なぜ Fama-Macbeth を 2 段階回帰モデルと呼ぶのかということです。 以下の手順を見てみましょう。
- このステップでは、時系列アプローチを使用して、1つまたは複数のリスク要因に対してすべての資産のリターンを回帰します。
ベータ、ファクター負荷、またはファクターエクスポージャーとして知られるすべてのファクターへのリターンエクスポージャーを取得します。 - このステップでは、クロスセクション アプローチを使用して、前のステップ (ステップ 1) で取得した資産
ベータに対してすべての資産のリターンを回帰します。 ここで、すべての要因のリスク プレミアムを取得します。
Fama と Macbeth によると、すべてのリスク要因に対するユニット エクスポージャーの経時的な予測プレミアムは、要素ごとに 1 回係数を平均することによって計算されます。
Python で Fama-Macbeth 回帰を実装する手順
2018年秋のFama-Macbethのライブラリ事情を反映したアップデートによると、Pythonのpandasモジュールからfama_macbethが削除されてしばらく経ちました。
では、Python を使用している場合、どのように Fama-Macbeth を実装できますか? このチュートリアルでは、次の 2つのアプローチを 1つずつ学習します。
- Python バージョン 3 を使用する場合、
LinearModelsでfama_macbeth関数を使用できます。 LinearModelsを使用する必要がない場合、または Python バージョン 2 を使用している場合、Fama-Macbeth の実装を作成するのがおそらく最良のケースです。
LinearModels を使用して Python で Fama-Macbeth 回帰を実装する
-
必要なモジュールとライブラリをインポートします。
import pandas as pd import numpy as np from pathlib import Path from linearmodels.asset_pricing import LinearFactorModel from statsmodels.api import OLS, add_constant import matplotlib.pyplot as plt import pandas_datareader.data as web import seaborn as sns import warnings warnings.filterwarnings("ignore") sns.set_style("whitegrid")まず、
LinearModelsを使用して Fama-Macbeth を実装するために必要なすべてのモジュールとライブラリをインポートします。 それらすべての簡単な説明を以下に示します。- データ フレームを操作するために
pandasをインポートし、配列のためにnumpyをインポートします。 pathlibは、この特定のスクリプトをPathオブジェクトに配置することにより、指定されたファイルへのパスを作成します。linearmodels.asset_pricingからLinearFactorModelをインポートします。 線形ファクター モデル は、線形方程式で記述される関係を使用して、資産の収益を限定/制限された数のファクターの値に関連付けます。- 次に、線形回帰モデルを評価するために
OLSをインポートし、配列に 1 の列を追加するためにadd_constantをインポートします。statsmodelsの詳細については こちら をご覧ください。 - その後、
pandas_datareaderをインポートして、pandasで使用する最新のリモート データにアクセスします。 さまざまなpandasバージョンで動作します。 - データのプロットと視覚化の目的で
matplotlibとseabornライブラリをインポートします。 Warningsをインポートして、警告を無視するfilterwarnings()メソッドを使用します。- 最後に、プロットの一般的なスタイルを制御するパラメータを設定する
seabornモジュールのset_style()メソッドを使用します。
- データ フレームを操作するために
-
リモートのリスク要因と研究ポートフォリオ データセットにアクセスします。
ff_research_data_five_factor = "F-F_Research_Data_5_Factors_2x3" ff_factor_dataset = web.DataReader( ff_research_data_five_factor, "famafrench", start="2010", end="2017-12" )[0]ここでは、ウェブサイト で入手できる更新されたリスク要因と研究ポートフォリオ データセット (Fama-French の 5つの要因) を使用して、上記のように
2010 年から 2017 年の月次頻度に戻しています。DataReader()を使用して、指定されたインターネット リソースから、コード フェンスのff_factor_datasetであるpandasデータ フレームにデータを抽出します。DataReader()はさまざまなソースをサポートしています。たとえば、Tiingo、IEX、Alpha Vantageなど、このページで読むことができます。print("OUTPUT for info(): \n") print(ff_factor_dataset.info())次に、
df.info()およびdf.describe()メソッドを使用します。ここで、info()は、列数、列のデータ型、列ラベル、範囲インデックス、メモリを含むデータ フレームの情報を出力します。 使用法、および各列のセル数 (null 以外の値)。以下の
info()によって生成された出力を確認できます。出力：
OUTPUT for info(): <class 'pandas.core.frame.DataFrame'> PeriodIndex: 96 entries, 2010-01 to 2017-12 Freq: M Data columns (total 6 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 Mkt-RF 96 non-null float64 1 SMB 96 non-null float64 2 HML 96 non-null float64 3 RMW 96 non-null float64 4 CMA 96 non-null float64 5 RF 96 non-null float64 dtypes: float64(6) memory usage: 5.2 KB None次に、次のように
describe()メソッドを使用します。print("OUTPUT for describe(): \n") print(ff_factor_dataset.describe())describe()メソッドは、データ フレームの統計概要を表示します。 この方法は Python シリーズにも使用できます。 この統計サマリーには、平均、中央値、度数、標準偏差、列のパーセンタイル値、および最小値と最大値が含まれます。以下の
describe()メソッドの出力を見つけることができます。出力：
OUTPUT for describe(): Mkt-RF SMB HML RMW CMA RF count 96.000000 96.000000 96.000000 96.000000 96.000000 96.000000 mean 1.158438 0.060000 -0.049271 0.129896 0.047708 0.012604 std 3.580012 2.300292 2.202912 1.581930 1.413033 0.022583 min -7.890000 -4.580000 -4.700000 -3.880000 -3.240000 0.000000 25% -0.917500 -1.670000 -1.665000 -1.075000 -0.952500 0.000000 50% 1.235000 0.200000 -0.275000 0.210000 0.010000 0.000000 75% 3.197500 1.582500 1.205000 1.235000 0.930000 0.010000 max 11.350000 7.040000 8.190000 3.480000 3.690000 0.090000 -
17 の業界ポートフォリオにアクセスし、リスク ファクター レートを差し引きます。
ff_industry_portfolio = "17_Industry_Portfolios" ff_industry_portfolio_dataset = web.DataReader( ff_industry_portfolio, "famafrench", start="2010", end="2017-12" )[0] ff_industry_portfolio_dataset = ff_industry_portfolio_dataset.sub( ff_factor_dataset.RF, axis=0 )ここでは、
DataReader()を使用して、毎月の頻度で 17 の業界ポートフォリオまたは資産にアクセスし、DataReader()によって返されたデータ フレームからリスクフリーレート (RF) を差し引きます。 なぜ？ これは、ファクター モデルが超過リターンで機能するためです。次に、
info()メソッドを使用して、チェーンされたデータ フレームに関する情報を取得します。print(ff_industry_portfolio_dataset.info())出力：
<class 'pandas.core.frame.DataFrame'> PeriodIndex: 96 entries, 2010-01 to 2017-12 Freq: M Data columns (total 17 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 Food 96 non-null float64 1 Mines 96 non-null float64 2 Oil 96 non-null float64 3 Clths 96 non-null float64 4 Durbl 96 non-null float64 5 Chems 96 non-null float64 6 Cnsum 96 non-null float64 7 Cnstr 96 non-null float64 8 Steel 96 non-null float64 9 FabPr 96 non-null float64 10 Machn 96 non-null float64 11 Cars 96 non-null float64 12 Trans 96 non-null float64 13 Utils 96 non-null float64 14 Rtail 96 non-null float64 15 Finan 96 non-null float64 16 Other 96 non-null float64 dtypes: float64(17) memory usage: 13.5 KB None同様に、
describe()メソッドを使用してこのデータ フレームを記述します。print(ff_industry_portfolio_dataset.describe())出力：
Food Mines Oil Clths Durbl Chems \ count 96.000000 96.000000 96.000000 96.000000 96.000000 96.000000 mean 1.046771 0.202917 0.597187 1.395833 1.151458 1.305000 std 2.800555 7.904401 5.480938 5.024408 5.163951 5.594161 min -5.170000 -24.380000 -11.680000 -10.000000 -13.160000 -17.390000 25% -0.785000 -5.840000 -3.117500 -1.865000 -2.100000 -1.445000 50% 0.920000 -0.435000 0.985000 1.160000 1.225000 1.435000 75% 3.187500 5.727500 4.152500 3.857500 4.160000 4.442500 max 6.670000 21.940000 15.940000 17.190000 16.610000 18.370000 Cnsum Cnstr Steel FabPr Machn Cars \ count 96.000000 96.000000 96.000000 96.000000 96.000000 96.000000 mean 1.186979 1.735521 0.559167 1.350521 1.217708 1.279479 std 3.142989 5.243314 7.389679 4.694408 4.798098 5.719351 min -7.150000 -14.160000 -20.490000 -11.960000 -9.070000 -11.650000 25% -0.855000 -2.410000 -4.395000 -1.447500 -2.062500 -1.245000 50% 1.465000 2.175000 0.660000 1.485000 1.525000 0.635000 75% 3.302500 5.557500 4.212500 3.837500 4.580000 4.802500 max 8.260000 15.510000 21.350000 17.660000 14.750000 20.860000 Trans Utils Rtail Finan Other count 96.000000 96.000000 96.000000 96.000000 96.000000 mean 1.463750 0.896458 1.233958 1.248646 1.290938 std 4.143005 3.233107 3.512518 4.839150 3.697608 min -8.560000 -6.990000 -9.180000 -11.140000 -7.890000 25% -0.810000 -0.737500 -0.952500 -1.462500 -1.090000 50% 1.480000 1.240000 0.865000 1.910000 1.660000 75% 4.242500 2.965000 3.370000 4.100000 3.485000 max 12.980000 7.840000 12.440000 13.410000 10.770000 -
超過リターンを計算します。
超過リターンの計算に進む前に、さらにいくつかの手順を実行する必要があります。
data_store = Path("./data/assets.h5") wiki_prices_df = pd.read_csv( "./dataset/wiki_prices.csv", parse_dates=["date"], index_col=["date", "ticker"], infer_datetime_format=True, ).sort_index() us_equities_data_df = pd.read_csv("./dataset/us_equities_data.csv") with pd.HDFStore(data_store) as hdf_store: hdf_store.put("quandl/wiki/prices", wiki_prices_df) with pd.HDFStore(data_store) as hdf_store: hdf_store.put("us_equities/stocks", us_equities_data_df.set_index("ticker"))まず、現在のディレクトリにある
dataフォルダーにassets.h5ファイルを作成します。 次に、read_csv()メソッドを使用して、同じディレクトリからdatasetフォルダー内のwiki_prices.csvおよびus_equities_data.csvファイルを読み取ります。その後、上記のように
HDFStore()を使用してHDF5形式でデータを保存します。with pd.HDFStore("./data/assets.h5") as hdf_store: prices = hdf_store["/quandl/wiki/prices"].adj_close.unstack().loc["2010":"2017"] equities = hdf_store["/us_equities/stocks"].drop_duplicates() sectors = equities.filter(prices.columns, axis=0).sector.to_dict() prices = prices.filter(sectors.keys()).dropna(how="all", axis=1) returns_df = prices.resample("M").last().pct_change().mul(100).to_period("M") returns_df = returns.dropna(how="all").dropna(axis=1) print(returns_df.info())上記のコード フェンスでは、
assets.h5ファイルに保存した/quandl/wiki/pricesと/us_equities/stocksを読み取り、それらをpricesおよびequities変数に保存しました。次に、
pricesとequitiesにいくつかのフィルターを適用し、リサンプリングして、欠落している値を削除します。 最後に、info()メソッドを使用して、データ フレームreturns_dfの情報を出力します。 以下の出力を確認できます。出力：
<class 'pandas.core.frame.DataFrame'> PeriodIndex: 95 entries, 2010-02 to 2017-12 Freq: M Columns: 1986 entries, A to ZUMZ dtypes: float64(1986) memory usage: 1.4 MB Noneここで、次のコードを実行してデータを整列させます。
ff_factor_dataset = ff_factor_dataset.loc[returns_df.index] ff_industry_portfolio_dataset = ff_industry_portfolio_dataset.loc[returns_df.index] print(ff_factor_dataset.describe())出力：
Mkt-RF SMB HML RMW CMA RF count 95.000000 95.000000 95.000000 95.000000 95.000000 95.000000 mean 1.206000 0.057053 -0.054316 0.144632 0.043368 0.012737 std 3.568382 2.312313 2.214041 1.583685 1.419886 0.022665 min -7.890000 -4.580000 -4.700000 -3.880000 -3.240000 0.000000 25% -0.565000 -1.680000 -1.670000 -0.880000 -0.965000 0.000000 50% 1.290000 0.160000 -0.310000 0.270000 0.010000 0.000000 75% 3.265000 1.605000 1.220000 1.240000 0.940000 0.010000 max 11.350000 7.040000 8.190000 3.480000 3.690000 0.090000これで、超過リターンを計算する準備が整いました。
excess_returns_df = returns_df.sub(ff_factor_dataset.RF, axis=0) excess_returns_df = excess_returns_df.clip( lower=np.percentile(excess_returns_df, 1), upper=np.percentile(excess_returns_df, 99), ) excess_returns_df.info()上記のコード ブロックでは、
ff_factor_datasetからリスク ファクターを減算し、返されたデータ フレームをexcess_returns_dfに保存します。次に、指定された入力しきい値で値をトリミングする
.clip()メソッドを使用します。 最後に、info()を使用して、excess_returns_dfデータ フレームの情報を出力します。出力：
<class 'pandas.core.frame.DataFrame'> PeriodIndex: 95 entries, 2010-02 to 2017-12 Freq: M Columns: 1986 entries, A to ZUMZ dtypes: float64(1986) memory usage: 1.4 MBFama-Macbeth 回帰の最初の段階に進む前に、
.drop()メソッドを使用します。このメソッドは、列または行から指定されたラベルを削除します。axis=1は列からのドロップを意味し、axis=0は行からのドロップを意味することに注意してください。ff_factor_dataset = ff_factor_dataset.drop("RF", axis=1) print(ff_factor_dataset.info())出力：
<class 'pandas.core.frame.DataFrame'> PeriodIndex: 95 entries, 2010-02 to 2017-12 Freq: M Data columns (total 5 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 Mkt-RF 95 non-null float64 1 SMB 95 non-null float64 2 HML 95 non-null float64 3 RMW 95 non-null float64 4 CMA 95 non-null float64 dtypes: float64(5) memory usage: 4.5 KB None -
Fama-Macbeth 回帰ステップ 1: ファクター エクスポージャーを実装します。
betas = [] for industry in ff_industry_portfolio_dataset: step_one = OLS( endog=ff_industry_portfolio_dataset.loc[ff_factor_dataset.index, industry], exog=add_constant(ff_factor_dataset), ).fit() betas.append(step_one.params.drop("const")) betas = pd.DataFrame( betas, columns=ff_factor_dataset.columns, index=ff_industry_portfolio_dataset.columns, ) print(betas.info())上記のコード スニペットは、Fama-Macbeth 回帰の最初のステップを実装し、17 因子の負荷推定値にアクセスします。 ここでは、
OLS()を使用して線形回帰モデルを評価し、add_constant()を使用して配列に 1 の列を追加します。出力：
<class 'pandas.core.frame.DataFrame'> Index: 17 entries, Food to Other Data columns (total 5 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 Mkt-RF 17 non-null float64 1 SMB 17 non-null float64 2 HML 17 non-null float64 3 RMW 17 non-null float64 4 CMA 17 non-null float64 dtypes: float64(5) memory usage: 1.3+ KB None -
Fama-Macbeth 回帰ステップ 2: リスクプレミアムを実装します。
lambdas = [] for period in ff_industry_portfolio_dataset.index: step_two = OLS( endog=ff_industry_portfolio_dataset.loc[period, betas.index], exog=betas ).fit() lambdas.append(step_two.params) lambdas = pd.DataFrame( lambdas, index=ff_industry_portfolio_dataset.index, columns=betas.columns.tolist() ) print(lambdas.info())2 番目のステップでは、ファクター負荷に関するポートフォリオのクロスセクションの
期間リターンの 96 回の回帰を実行します。出力：
<class 'pandas.core.frame.DataFrame'> PeriodIndex: 95 entries, 2010-02 to 2017-12 Freq: M Data columns (total 5 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 Mkt-RF 95 non-null float64 1 SMB 95 non-null float64 2 HML 95 non-null float64 3 RMW 95 non-null float64 4 CMA 95 non-null float64 dtypes: float64(5) memory usage: 6.5 KB None結果は次のように視覚化できます。
window = 24 # here 24 is the number of months axis1 = plt.subplot2grid((1, 3), (0, 0)) axis2 = plt.subplot2grid((1, 3), (0, 1), colspan=2) lambdas.mean().sort_values().plot.barh(ax=axis1) lambdas.rolling(window).mean().dropna().plot( lw=1, figsize=(14, 5), sharey=True, ax=axis2 ) sns.despine() plt.tight_layout()出力：
-
LinearModelsモジュールを使用して Fama-Macbeth 回帰を実装します。model = LinearFactorModel( portfolios=ff_industry_portfolio_dataset, factors=ff_factor_dataset ) result = model.fit() print(result)ここでは、
LinearModelsを使用して、次の出力を生成する 2 ステップの Fama-Macbeth 手順を実装します。print(result)の代わりにprint(result.full_summary)を使用して完全な要約を得ることができます。出力：
LinearFactorModel Estimation Summary ================================================================================ No. Test Portfolios: 17 R-squared: 0.6879 No. Factors: 5 J-statistic: 15.619 No. Observations: 95 P-value 0.2093 Date: Mon, Oct 24 2022 Distribution: chi2(12) Time: 20:53:52 Cov. Estimator: robust Risk Premia Estimates ============================================================================== Parameter Std. Err. T-stat P-value Lower CI Upper CI ------------------------------------------------------------------------------ Mkt-RF 1.2355 0.4098 3.0152 0.0026 0.4324 2.0386 SMB 0.0214 0.8687 0.0246 0.9804 -1.6813 1.7240 HML -1.1140 0.6213 -1.7931 0.0730 -2.3317 0.1037 RMW -0.2768 0.8133 -0.3403 0.7336 -1.8708 1.3172 CMA -0.5078 0.5666 -0.8962 0.3701 -1.6183 0.6027 ============================================================================== Covariance estimator: HeteroskedasticCovariance See full_summary for complete results
Python で Fama-Macbeth 回帰を実装するための代替アプローチ
LinearModels または Python バージョン 2 を使用したくない場合は、このアプローチを使用できます。
-
モジュールとライブラリをインポートします。
import pandas as pd import numpy as np import statsmodels.formula.api as smpandasをインポートしてデータ フレームを操作し、numpyを配列で操作し、statsmodels.formula.apiをインポートします。statsmodels.formula.apiは、数式文字列とデータ フレームを介してモデルを指定するための便利なインターフェイスです。 -
データセットの読み取りとクエリ。
以下に示すように、パネルにFama-French業界の資産/ポートフォリオがあると仮定します（たとえば、過去の
ベータとリターンなどのいくつかの変数を計算して、x変数として使用します）data_df = pd.read_csv("industry_data.csv", parse_dates=["caldt"]) data_df.query("caldt == '1995-07-01'")出力：
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 -
groupby()を使用して、クロスセクション回帰モデルを月ごとに計算します。def ols_coefficient(x, formula): return sm.ols(formula, data=x).fit().params gamma_df = data_df.groupby("caldt").apply( ols_coefficient, "ret ~ 1 + beta + r12to2 + r36to13" ) gamma_df.head()ここで、
groupbyを使用するのは、Fama-Macbeth が正確な横断回帰モデルを月ごとに計算する必要があるためです。 データ フレーム (groupbyによって返される) とpasty式を受け入れる関数を作成できます。 次に、モデルに適合し、パラメーターの推定値を返します。出力：
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 -
平均値と平均値の標準誤差を計算します。
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"]] fm_summary(gamma_df)次に、平均値、t 検定 (任意の統計を使用できます)、および平均値の標準誤差を計算します。これは上記のようなものです。
出力：
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 -
速度を改善し、
fama_macbeth関数を使用します。def ols_np(dataset, y_var, x_var): gamma_df, _, _, _ = np.linalg.lstsq(dataset[x_var], dataset[y_var], rcond=None) return pd.Series(gamma_df)効率を重視する場合、このステップは重要です。 もしそうなら、
statsmodelsからnumpy.linalg.lstsqに切り替えることができます。ols推定を実行するには、上記のような関数を書くことができます。 これらの行列のランクをチェックするために何もしていないことに注意してください。古い
pandasバージョンを使用している場合は、次のように動作するはずです。pandasでfama_macbethを使用した例を見てみましょう。print(data_df) fm = pd.fama_macbeth(y=data_df["y"], x=data_df[["x"]]) print(fm)ここで、以下の構造を観察してください。
fama_macbethは、x-varとy-varに、インデックスの最初の変数としてdateと 2 番目の変数としてstock/firm/entityid を持つマルチインデックスを持たせたいと考えています。出力：
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 ----------------------Summary of the Fama-Macbeth Analysis------------------------- Formula: Y ~ x + intercept # betas : 3 ----------------------Summary of the 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 the Summary---------------------------------上記のコード フェンスで行ったように、
fmを出力するだけでfm.summaryを呼び出していることに注意してください。 さらに、fama_macbethは自動的にinterceptを追加し、x-varはデータ フレームでなければなりません。interceptを使用したくない場合は、次のようにできます。fm = pd.fama_macbeth(y=data_df["y"], x=data_df[["x"]], intercept=False)
