linearmodels#

in English or the language of your choice.

import pandas as pd
import py4macro
import wooldridge
from linearmodels.panel.data import PanelData
from linearmodels.panel import FirstDifferenceOLS

# 警告メッセージを非表示
import warnings
warnings.filterwarnings("ignore")

説明#

linearmodelsstatsmodelsを補完する目的として開発されている。主に,パネルデータ,操作変数法を使った推定法やGMMを扱う場合には非常に重宝するパッケージである。しかし,linearmodelsstatsmodelsの両方を使う上で以下の点に注意する必要がある。

  • 推定結果などのメソッドや属性が共通化されているわけではない。次の表に3つの例を挙げる。

推定結果の
表を表示

残差を
表示する
メソッド

標準誤差を
取得する
属性

statsmodels

.summary()

.resid

.bse

linearmodels

.summary

.resids

.std_errors

  • statsmodelslinearmodelsも回帰式を文字列で指定できるが,定数項を指定する方法が異なる。

    • statsmodelsでは,定数項は自動的に追加される,定数項を省く場合は-1を追加する。

    • linearmodelsでは,定数項は自動的に追加されない。定数項を入れる場合は1を追加する。

  • fit()メソッドの挙動も共通化されていない。

    • linearmodelsfit()に何のオプションも指定せずにOLS推定すると係数の推定量は同じだが,標準誤差や\(t\)値などが異なる。同じにするためには次のように2つのオプションを設定しなくてはならない。

    .fit(cov_type='unadjusted', debiased=True)
    
    • cov_typeは不均一分散頑健共分散行列推定のオプション

      • デフォルトはrobust(不均一分散頑健的共分散行列推定)でstatsmodelsHC1と等しい。

    • debiasedは共分散行列推定の自由度のオプション(小標本の場合の調整)

      • デフォルトはFalse

(注意)

以下では.fit()のオプションは指定せず,デフォルトのまま議論を続ける。

Note

linearmodelsを事前にインストールする必要があるが,次のコマンドでインストールすることができる。

pip install linearmodels

以下ではlinearmodelsを使うが,そのためにはDataFrameMultiIndexに変換する必要がある。以下では,まずMultiIndexについて説明し,その後にlinearmodelsにあるPanelData関数について説明する

PandasMultiIndex#

説明#

パネル・データを扱うために必要なPandasMultiIndexについて説明する。MultiIndexとは行や列のラベルが階層的になったDataFrameSeriesを指す。以下では,DataFrameの行におけるMultiIndexを説明する。

まずデータを読み込む。

# url の設定
url = 'https://raw.githubusercontent.com/Haruyama-KobeU/Haruyama-KobeU.github.io/master/data/data4.csv'

# 読み込み
df = pd.read_csv(url)
df
year country gdp inv con pop
0 2000 Australia 90 30.0 50.0 10
1 2001 Australia 100 40.0 80.0 11
2 2002 Australia 120 NaN 70.0 16
3 2000 Japan 100 20.0 80.0 8
4 2001 Japan 95 25.0 70.0 9
5 2002 Japan 93 21.0 72.0 10
6 2000 UK 100 30.0 70.0 11
7 2001 UK 110 39.0 71.0 12
8 2002 UK 115 55.0 NaN 14

行・列ともにラベルの階層は1つずつとなっている。set_index()を使い行にMultiIndexを作成するが,引数に

\[\left[\text{第0id},\text{第1id}\right]\]

とし階層インデックス化する。ここでパネル・データ分析をする上で以下のルールに従うことにする。

  • 第0id:観察単位(例えば,消費者,企業,国)

  • 第1id:時間(例えば,年,四半期)

次の例ではcountryyearの行をそれぞれ第0インデックス,第1インデックスに指定する。

df = df.set_index(['country', 'year'])#.sort_index()
df
gdp inv con pop
country year
Australia 2000 90 30.0 50.0 10
2001 100 40.0 80.0 11
2002 120 NaN 70.0 16
Japan 2000 100 20.0 80.0 8
2001 95 25.0 70.0 9
2002 93 21.0 72.0 10
UK 2000 100 30.0 70.0 11
2001 110 39.0 71.0 12
2002 115 55.0 NaN 14

階層インデックスが綺麗に並んでいるが,元のデータの並び方によっては階層インデックスが期待通りに並ばない場合がありえる。その場合は,メソッドsort_index()を使うと良いだろう。


MultiIndexを解除するにはメソッド.reset_index()を使う。

df.reset_index()
country year gdp inv con pop
0 Australia 2000 90 30.0 50.0 10
1 Australia 2001 100 40.0 80.0 11
2 Australia 2002 120 NaN 70.0 16
3 Japan 2000 100 20.0 80.0 8
4 Japan 2001 95 25.0 70.0 9
5 Japan 2002 93 21.0 72.0 10
6 UK 2000 100 30.0 70.0 11
7 UK 2001 110 39.0 71.0 12
8 UK 2002 115 55.0 NaN 14

要素,行,列の抽出#

MultiIndexのまま要素・列・行の抽出およびスライシングには様々な方法があり,複雑である。特に,スライシングをしたい場合,一番簡単なのはreset_index()で通常のDataFrameに戻し,スライシングし新たなDataFrameを作成するだけでも十分であろう。

以下では,.loc[]を使いMultiIndexのままでの抽出方法について簡単に説明する。その際,以下のルールは変わらない。

\[.\text{loc}\left[\text{行の指定},\text{列の指定}\right]\]

ただ,行の指定にリストやタプルを使うことになる(列の指定も同じ)。

他の方法についてはこのサイトこのサイトが参考になる。

1つの観察単位の抽出#

1つの要素を抽出する場合はタプルを使う。例えば,日本の2001年のgdpを抽出したい場合。

df.loc[('Japan',2001), 'gdp']
95

行の抽出#

上の例で列の指定を:にすると,指定した行に対して全ての列を抽出できる。

df.loc[('Japan',2001), :]
gdp    95.0
inv    25.0
con    70.0
pop     9.0
Name: (Japan, 2001), dtype: float64

この場合,特定の列に対してスライシングも可能。

df.loc[('Japan',2001), 'gdp':'con']
gdp    95.0
inv    25.0
con    70.0
Name: (Japan, 2001), dtype: float64

指定した行に対して個別に複数列を抽出したい場合は,タプルを使う。

df.loc[('Japan',2001), ('gdp','con')]
gdp    95.0
con    70.0
Name: (Japan, 2001), dtype: float64

複数行の抽出にはリストで指定する。

df.loc[(['Japan','UK'],[2001,2002]), :]
gdp inv con pop
country year
Japan 2001 95 25.0 70.0 9
2002 93 21.0 72.0 10
UK 2001 110 39.0 71.0 12
2002 115 55.0 NaN 14

第0インデックスの観察単位の全て#

第0インデックスにある,ある観察単位の全てのデータだけを抽出したい場合は,通常のPandasの場合と同じ。

df.loc['Japan', :]
gdp inv con pop
year
2000 100 20.0 80.0 8
2001 95 25.0 70.0 9
2002 93 21.0 72.0 10

複数の場合。

df.loc[['Japan','UK'], :]
gdp inv con pop
country year
Japan 2000 100 20.0 80.0 8
2001 95 25.0 70.0 9
2002 93 21.0 72.0 10
UK 2000 100 30.0 70.0 11
2001 110 39.0 71.0 12
2002 115 55.0 NaN 14

列の抽出#

通常のPandasと同じ。Seriesを返す場合。

df.loc[:,'gdp']
country    year
Australia  2000     90
           2001    100
           2002    120
Japan      2000    100
           2001     95
           2002     93
UK         2000    100
           2001    110
           2002    115
Name: gdp, dtype: int64

[]を使うと,DataFrameとして抽出できる。

df.loc[:,['gdp']]
gdp
country year
Australia 2000 90
2001 100
2002 120
Japan 2000 100
2001 95
2002 93
UK 2000 100
2001 110
2002 115

複数列抽出の場合。

df.loc[:,['gdp','inv']]
gdp inv
country year
Australia 2000 90 30.0
2001 100 40.0
2002 120 NaN
Japan 2000 100 20.0
2001 95 25.0
2002 93 21.0
UK 2000 100 30.0
2001 110 39.0
2002 115 55.0

スライシングも使える。

df.loc[:,'gdp':'con']
gdp inv con
country year
Australia 2000 90 30.0 50.0
2001 100 40.0 80.0
2002 120 NaN 70.0
Japan 2000 100 20.0 80.0
2001 95 25.0 70.0
2002 93 21.0 72.0
UK 2000 100 30.0 70.0
2001 110 39.0 71.0
2002 115 55.0 NaN

第1インデックスのある年だけの抽出#

一番簡単な方法はreset_index()を使い今まで習った関数を使う。

df.reset_index().query('year == 2000')
country year gdp inv con pop
0 Australia 2000 90 30.0 50.0 10
3 Japan 2000 100 20.0 80.0 8
6 UK 2000 100 30.0 70.0 11

複数年の場合。

df.reset_index().query('year in [2000,2002]')
country year gdp inv con pop
0 Australia 2000 90 30.0 50.0 10
2 Australia 2002 120 NaN 70.0 16
3 Japan 2000 100 20.0 80.0 8
5 Japan 2002 93 21.0 72.0 10
6 UK 2000 100 30.0 70.0 11
8 UK 2002 115 55.0 NaN 14

上と同じ結果。

df.reset_index().query('year not in [2001]')
country year gdp inv con pop
0 Australia 2000 90 30.0 50.0 10
2 Australia 2002 120 NaN 70.0 16
3 Japan 2000 100 20.0 80.0 8
5 Japan 2002 93 21.0 72.0 10
6 UK 2000 100 30.0 70.0 11
8 UK 2002 115 55.0 NaN 14

linearmodelsPanelData#

linearmodelsではMultiIndex化されたDataFrameをそのまま読み込み推定することができる。一方で,linearmodelsの関数PanelDataを使いMultiIndex化されたDataFramePanelDataオブジェクトに変換すると分析に必要な計算を簡単にできるようになる。必須ではないが,知っておいて損はしない関数である。

まずdfPanelDataオブジェクトに変換する。

dfp = PanelData(df)
dfp
PanelData
gdp inv con pop
country year
Australia 2000 90.0 30.0 50.0 10.0
2001 100.0 40.0 80.0 11.0
2002 120.0 NaN 70.0 16.0
Japan 2000 100.0 20.0 80.0 8.0
2001 95.0 25.0 70.0 9.0
2002 93.0 21.0 72.0 10.0
UK 2000 100.0 30.0 70.0 11.0
2001 110.0 39.0 71.0 12.0
2002 115.0 55.0 NaN 14.0

属性とメソッド#

まず,py4macroモジュールに含まれるsee()関数を使い,dfpの属性とメソッドに何があるかを確認する。

py4macro.see(dfp)
.copy               .count              .dataframe          .demean
.drop               .dummies            .entities           .entity_ids
.first_difference   .general_demean     .index              .isnull
.mean               .ndim               .nentity            .nobs
.nvar               .panel              .shape              .time
.time_ids           .values2d           .values3d           .vars

主なものについて説明する。

属性shapeは,PanelDataの変数の数を表示する。以下が返り値の内容である。

(変数の数, 時間の観測値の数, 観察単位の数)
dfp.shape
(4, 3, 3)
  • 変数の数:4(列にある変数)

  • 時間の観測値の数:3(年)

  • 観察単位の数:3(国)

メソッド.mean()を使うと、変数の観察単位毎の平均のDataFrameが返される。

dfp.mean()
gdp inv con pop
country
Australia 103.333333 35.000000 66.666667 12.333333
Japan 96.000000 22.000000 74.000000 9.000000
UK 108.333333 41.333333 70.500000 12.333333

引数にtimeを指定すると、変数の時間毎の平均が返される。

dfp.mean('time')
gdp inv con pop
year
2000 96.666667 26.666667 66.666667 9.666667
2001 101.666667 34.666667 73.666667 10.666667
2002 109.333333 38.000000 71.000000 13.333333

メソッドdemean()は、変数の平均からの乖離が返される。即ち、変数\(x\)の平均が\(\bar{x}\)とすると、\(x-\bar{x}\)が返される。

dfp.demean()
PanelData
gdp inv con pop
country year
Australia 2000 -13.333333 -5.000000 -16.666667 -2.333333
2001 -3.333333 5.000000 13.333333 -1.333333
2002 16.666667 NaN 3.333333 3.666667
Japan 2000 4.000000 -2.000000 6.000000 -1.000000
2001 -1.000000 3.000000 -4.000000 0.000000
2002 -3.000000 -1.000000 -2.000000 1.000000
UK 2000 -8.333333 -11.333333 -0.500000 -1.333333
2001 1.666667 -2.333333 0.500000 -0.333333
2002 6.666667 13.666667 NaN 1.666667

first_difference()は変数の1階差分(\(x_t-x_{t-1}\))が返される。

dfp.first_difference()
PanelData
gdp inv con pop
country year
Australia 2001 10.0 10.0 30.0 1.0
Japan 2001 -5.0 5.0 -10.0 1.0
2002 -2.0 -4.0 2.0 1.0
UK 2001 10.0 9.0 1.0 1.0

上の例ではNaNがあるためAustraliaUKの行は1つしかない。


(注意)

DataFrameのメソッドはPanelDataオブジェクトには使えない。

従って,DataFrameのメソッド(例えば,行や列の抽出)を使う場合,DataFrameに変換する必要がある。その場合,PanelDataオブジェクトの属性.dataframeを使うことができる。

dfp.dataframe.loc['Japan',:]
gdp inv con pop
year
2000 100.0 20.0 80.0 8.0
2001 95.0 25.0 70.0 9.0
2002 93.0 21.0 72.0 10.0

Balanced/Unbalancedの確認#

データセットには欠損値がある場合がある。観察単位数が\(N\)で時間の観測値の数が\(T\)の場合,観測値の数は\(n=N\times T\)となるが,次の2つを区別する。

  • balanced panel data:\(n=N\times T\)(観察単位に対して全ての期間の全ての変数に欠損値がない)

  • unbalanced panel data:\(n<N\times T\)(欠損値がある)

balanced か unbalancedかは以下のコードで確認できる。まず,属性isnullを使う。

dfp.isnull
country    year
Australia  2000    False
           2001    False
           2002     True
Japan      2000    False
           2001    False
           2002    False
UK         2000    False
           2001    False
           2002     True
dtype: bool

それぞれの行にNaNがあればTrueを、なければFalseを返す。次にTrue/Falseを逆転させるために~を使う。

~dfp.isnull
country    year
Australia  2000     True
           2001     True
           2002    False
Japan      2000     True
           2001     True
           2002     True
UK         2000     True
           2001     True
           2002    False
dtype: bool

Trueの行にはNaNはなく、Falseの行にNaNがある。行数が多い場合はメソッドall()が便利である。all()は列に対して全ての要素がTrueの場合のみTrueを返す。

(~dfp.isnull).all()
False

Falseなので unbalanced panel data ということが確認できた。

1階差分推定(再考)#

ここではlinearmodelsを使い,以前行った1階差分推定を再考する。データcrime4を使う。

crime4 = wooldridge.data('crime4')
crime4.head()
county year crmrte prbarr prbconv prbpris avgsen polpc density taxpc ... lpctymle lpctmin clcrmrte clprbarr clprbcon clprbpri clavgsen clpolpc cltaxpc clmix
0 1 81 0.039885 0.289696 0.402062 0.472222 5.61 0.001787 2.307159 25.697630 ... -2.433870 3.006608 NaN NaN NaN NaN NaN NaN NaN NaN
1 1 82 0.038345 0.338111 0.433005 0.506993 5.59 0.001767 2.330254 24.874252 ... -2.449038 3.006608 -0.039376 0.154542 0.074143 0.071048 -0.003571 -0.011364 -0.032565 0.030857
2 1 83 0.030305 0.330449 0.525703 0.479705 5.80 0.001836 2.341801 26.451443 ... -2.464036 3.006608 -0.235316 -0.022922 0.193987 -0.055326 0.036879 0.038413 0.061477 -0.244732
3 1 84 0.034726 0.362525 0.604706 0.520104 6.89 0.001886 2.346420 26.842348 ... -2.478925 3.006608 0.136180 0.092641 0.140006 0.080857 0.172213 0.026930 0.014670 -0.027331
4 1 85 0.036573 0.325395 0.578723 0.497059 6.55 0.001924 2.364896 28.140337 ... -2.497306 3.006608 0.051825 -0.108054 -0.043918 -0.045320 -0.050606 0.020199 0.047223 0.172125

5 rows × 59 columns

countyyearを使いMultiIndex化する。

crime4 = crime4.set_index(['county','year'])
crime4.head()
crmrte prbarr prbconv prbpris avgsen polpc density taxpc west central ... lpctymle lpctmin clcrmrte clprbarr clprbcon clprbpri clavgsen clpolpc cltaxpc clmix
county year
1 81 0.039885 0.289696 0.402062 0.472222 5.61 0.001787 2.307159 25.697630 0 1 ... -2.433870 3.006608 NaN NaN NaN NaN NaN NaN NaN NaN
82 0.038345 0.338111 0.433005 0.506993 5.59 0.001767 2.330254 24.874252 0 1 ... -2.449038 3.006608 -0.039376 0.154542 0.074143 0.071048 -0.003571 -0.011364 -0.032565 0.030857
83 0.030305 0.330449 0.525703 0.479705 5.80 0.001836 2.341801 26.451443 0 1 ... -2.464036 3.006608 -0.235316 -0.022922 0.193987 -0.055326 0.036879 0.038413 0.061477 -0.244732
84 0.034726 0.362525 0.604706 0.520104 6.89 0.001886 2.346420 26.842348 0 1 ... -2.478925 3.006608 0.136180 0.092641 0.140006 0.080857 0.172213 0.026930 0.014670 -0.027331
85 0.036573 0.325395 0.578723 0.497059 6.55 0.001924 2.364896 28.140337 0 1 ... -2.497306 3.006608 0.051825 -0.108054 -0.043918 -0.045320 -0.050606 0.020199 0.047223 0.172125

5 rows × 57 columns

次にPanelDataオブジェクトに変換しデータの特徴を調べる。

crime4p = PanelData(crime4)
crime4p.shape
(57, 7, 90)
  • 57: 変数の数

  • 7: 時間の観測値の数(年次データなので7年間)

  • 90:観察単位の数(人数)

次にbalanced もしくは unbalanced data set かを確認する。

(~crime4p.isnull).all()
False

Unbalancedのデータセットだと確認できた。


実際に回帰式を書くことにする。使い方はstatsmodelsと似ている。

  • FirstDifferenceOLSモジュールの関数.from_formulaを使い次のように引数を指定する。

\[\text{.from_formula}(\text{回帰式}, \text{データ})\]
  • 定数項を入れることはできない仕様となっている。

  • ここでは,以前の推定結果と比べるために,ダミー変数d82を追加する。

formula = 'lcrmrte ~ d82 + d83 + d84 + d85 + d86 + d87 + lprbarr + \
                lprbconv + lprbpris + lavgsen + lpolpc'

1階差分モデルの設定(インスタンスの作成)

mod_dif = FirstDifferenceOLS.from_formula(formula, data=crime4)

statsmodelsと同じように,そこから得た結果にメソッド.fit()を使い計算し結果が返される。

res_dif = mod_dif.fit()

<結果の表示方法>

  1. res_difもしくはprint(res_dif)を実行。

  2. res_difには属性summaryが用意されているが、表示方法1と同じ内容が表示される。

  3. summaryには属性tablesがあり,2つの表がリストとして格納されている。

    • tables[0]:検定統計量の表(print()を使うと見やすくなる)

    • tables[1]:係数の推定値やp値などの表(print()を使うと見やすくなる)

print(res_dif.summary.tables[1])
                             Parameter Estimates                              
==============================================================================
            Parameter  Std. Err.     T-stat    P-value    Lower CI    Upper CI
------------------------------------------------------------------------------
d82            0.0077     0.0171     0.4522     0.6513     -0.0258      0.0412
d83           -0.0844     0.0235    -3.5998     0.0003     -0.1305     -0.0384
d84           -0.1247     0.0287    -4.3366     0.0000     -0.1811     -0.0682
d85           -0.1216     0.0331    -3.6670     0.0003     -0.1867     -0.0564
d86           -0.0863     0.0367    -2.3539     0.0189     -0.1584     -0.0143
d87           -0.0378     0.0400    -0.9455     0.3448     -0.1163      0.0407
lprbarr       -0.3275     0.0300    -10.924     0.0000     -0.3864     -0.2686
lprbconv      -0.2381     0.0182    -13.058     0.0000     -0.2739     -0.2023
lprbpris      -0.1650     0.0260    -6.3555     0.0000     -0.2161     -0.1140
lavgsen       -0.0218     0.0221    -0.9850     0.3250     -0.0652      0.0216
lpolpc         0.3984     0.0269     14.821     0.0000      0.3456      0.4512
==============================================================================

推定結果は以前のものと同じである。