固定効果モデル

パネルデータ¶

異なる時点で同じ個体についての情報を観測したデータをパネルデータ （panel data）という。

被説明変数 $Y$ と説明変数 $X = (X_1, \dots, X_K)$ に時点 $t$ と個体 $i$ の添字をつけて

\{ (Y_{it}, X_{it}) \}, \hspace{2em} i = 1, \dots, N, \hspace{1em} t = 1, \dots, T

(1)

のような形で表記される

import statsmodels.api as sm
import statsmodels.formula.api as smf

data = sm.datasets.get_rdataset("Grunfeld", package="plm").data
data

# Grunfeldは1935～1954年にかけてのアメリカの10の企業のbalanced panelデータ
# firm: 企業ID
# inv: 投資総額
# value: 企業価値
# capital: 資本ストック

Loading...

Pooled OLS¶

パネルデータ分析において、個体や時間による固定効果（fixed effect）を特に考慮しないで（固定効果が無いと仮定して）通常の重回帰モデルを用いたモデルをPooled OLSと呼ぶ。

Y_{it} = \boldsymbol{\beta} X_{it} + \varepsilon_{it}

(2)

pooled_ols = sm.OLS.from_formula(formula="inv ~ value + capital", data=data).fit()
pooled_ols.summary()

Loading...

固定効果モデル¶

パネルデータを分析する際に個体差による効果（個体固定効果）や時点ごとの固有の効果（時間固定効果）の影響を除くように作ったモデルのことを固定効果モデル（fixed effect model）という。

Y_{it} = \boldsymbol{\beta} X_{it} +\theta_i + \pi_t + \varepsilon_{it}

(3)

個体内の変動を使うため、固定効果推定量は「群内（within）推定量」と呼ばれたりもする。

one-way fixed effect model¶

時点や個体など、1つの固定効果に対処するモデル。一元配置固定効果モデル（one-way fixed effect model）などと呼ばれる。その推定量は「群内（within）」推定量や「固定効果（fixed effect）」推定量などとも呼ばれる。

個体固定効果モデル¶

以下のような個体固定効果モデルを考える。（単純化のため説明変数 $X$ は1つのみとする）

Y_{it} = \beta_0 + \beta_1 X_{it} +\theta_i + \varepsilon_{it} \hspace{2em} (i = 1, \dots, N, \ t=1,\dots, T)

(4)

パネルデータを用いることができる場合、以下の3つの方法によって個体固定効果（entity fixed effects） $\theta_i$ を除去することができる。

(1) "一回の階差モデル（first difference model）"によるOLS推定¶

(Y_{i,t+1} - Y_{it}) = (\beta_0 - \beta_0) + \beta_1 (X_{i,t+1}-X_{it}) + (\varepsilon_{i,t+1} - \varepsilon_{it})

(5)

記号を置き換えて、

\Delta Y_{it} = \beta_1 \Delta X_{it}+ \Delta \varepsilon_{it}

(6)

推定方法：
1. 説明変数、被説明変数それぞれ $t+1$ 期から $t$ 期を引く
2. 上の式をOLS推定する

import pandas as pd
data = sm.datasets.get_rdataset("Grunfeld", package="plm").data

deltas = []
data = data.sort_values(["firm", "year"])
for firm in data["firm"].unique():
    d = data.query(f"firm == {firm}").copy()
    delta = d - d.shift(1)
    delta["year"] = d["year"]
    delta["firm"] = firm
    deltas.append(delta)
delta = pd.concat(deltas).dropna().sort_values("firm").reset_index(drop=True)
delta

Loading...

first_diff = sm.OLS.from_formula(formula="inv ~ -1 + value + capital", data=delta).fit()
first_diff.summary()

Loading...

(2) “ $N-1$ 個のダミー説明変数”を用いたOLS推定¶

最小二乗ダミー変数推定（Least Squares Dummy Variables (LSDV) 推定）とも呼ばれる。

Y_{it} = \beta_0 + \beta_1 X_{it} + \gamma_2 D^{(2)}_i + \cdots + \gamma_N D^{(N)}_i + \varepsilon_{it} \\ \text{where } D^{(2)}_i = \begin{cases} 1 & \text{for } i = 2\\ 0 & \text{otherwise} \end{cases} \text{, etc.}

(7)

推定方法：
1. 個体ダミー変数（個体 $i$ に該当する場合に1、それ以外は0となるダミー変数） $D^{(2)}_i, \cdots, D^{(N)}_i$ を作成する
2. 上の式をOLS推定する

data = sm.datasets.get_rdataset("Grunfeld", package="plm").data

lsdv = sm.OLS.from_formula(
    formula="inv ~ value + capital + firm",
    data=data.assign(firm = data["firm"].astype("category")) # category型にすれば自動でダミー変数にしてくれる
).fit()
lsdv.summary()

Loading...

(3) ”平均差分法（Entity-demeaned）”を用いたOLS推定¶

\begin{align} \tilde{Y}_{it} &= \beta_1 \tilde{X}_{it} + \tilde{\varepsilon}_{it}, \\ \text{where } \tilde{Y}_{it} &= Y_{it} - \bar{Y}_i, \hspace{1em} \bar{Y}_i = \frac{1}{T} \sum^T_{t=1} Y_{it}\\ \tilde{X}_{it} &= X_{it} - \bar{X}_i, \hspace{1em} \bar{X}_i = \frac{1}{T} \sum^T_{t=1} X_{it}\\ \tilde{\varepsilon}_{it} &= \varepsilon_{it}- \bar{\varepsilon}_i, \hspace{1em} \bar{\varepsilon}_i = \frac{1}{T} \sum_{t=1}^T \varepsilon_{it} \end{align}

(8)

推定方法：
1. 説明変数・被説明変数について、変数から期間平均を引く
2. 上の式をOLS推定する
$N-1$ 個の個体ダミー説明変数による推定と同じ推定値が得られる
統計ソフトでは通常は平均差分法による推定が行われる

data = sm.datasets.get_rdataset("Grunfeld", package="plm").data
group = "firm"

rows = []
for _, d in data.groupby(group):
    for col in ["value", "inv", "capital"]:
        d[col] = (d[col] - d[col].mean())
    rows.append(d)
df = pd.concat(rows)

entity_demeaned = sm.OLS.from_formula(formula="inv ~ -1 + value + capital", data=df).fit()
entity_demeaned.summary()

Loading...

md = smf.mixedlm("inv ~ value + capital", data, groups=data["year"])
mdf = md.fit()
print(mdf.summary())

/usr/local/lib/python3.9/site-packages/statsmodels/base/model.py:607: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  warnings.warn("Maximum Likelihood optimization failed to "
/usr/local/lib/python3.9/site-packages/statsmodels/regression/mixed_linear_model.py:2201: ConvergenceWarning: Retrying MixedLM optimization with lbfgs
  warnings.warn(
/usr/local/lib/python3.9/site-packages/statsmodels/base/model.py:607: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  warnings.warn("Maximum Likelihood optimization failed to "
/usr/local/lib/python3.9/site-packages/statsmodels/regression/mixed_linear_model.py:2201: ConvergenceWarning: Retrying MixedLM optimization with cg
  warnings.warn(

         Mixed Linear Model Regression Results
========================================================
Model:            MixedLM Dependent Variable: inv       
No. Observations: 200     Method:             REML      
No. Groups:       20      Scale:              8912.8075 
Min. group size:  10      Log-Likelihood:     -1196.1133
Max. group size:  10      Converged:          No        
Mean group size:  10.0                                  
--------------------------------------------------------
            Coef.  Std.Err.   z    P>|z|  [0.025  0.975]
--------------------------------------------------------
Intercept  -42.714    9.501 -4.496 0.000 -61.335 -24.092
value        0.116    0.006 19.854 0.000   0.104   0.127
capital      0.231    0.025  9.136 0.000   0.181   0.280
Group Var    0.585                                      
========================================================

/usr/local/lib/python3.9/site-packages/statsmodels/base/model.py:607: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  warnings.warn("Maximum Likelihood optimization failed to "
/usr/local/lib/python3.9/site-packages/statsmodels/regression/mixed_linear_model.py:2207: ConvergenceWarning: MixedLM optimization failed, trying a different optimizer may help.
  warnings.warn(msg, ConvergenceWarning)
/usr/local/lib/python3.9/site-packages/statsmodels/regression/mixed_linear_model.py:2219: ConvergenceWarning: Gradient optimization failed, |grad| = 1.131047
  warnings.warn(msg, ConvergenceWarning)
/usr/local/lib/python3.9/site-packages/statsmodels/regression/mixed_linear_model.py:2262: ConvergenceWarning: The Hessian matrix at the estimated parameter values is not positive definite.
  warnings.warn(msg, ConvergenceWarning)

導出¶

固定効果のある

Y_{it} = \beta_0 + \beta_1 X_{it} +\theta_i + \varepsilon_{it}

(9)

というモデルを考えるとき、その推定量は

\newcommand{\argmin} attempting to redefine \argmin; use \renewcommand

(\hat{\beta}_0, \hat{\beta}_1, \hat{\theta}_1, \hat{\theta}_2, \dots, \hat{\theta}_N)
\newcommand{\argmin}{\mathop{\rm arg~min}\limits}
= \argmin_{\beta_0, \beta_1, \theta_1, \dots, \theta_N}
\sum^N_{i=1} \sum^T_{t=1} (Y_{it} - \beta_0 - \beta_1 X_{it} - \theta_i)^2

より、この最小化問題の一階条件は

\sum^N_{i=1} \sum^T_{t=1} X_{it} (Y_{it} - \hat{\beta}_0 - \hat{\beta}_1 X_{it} - \hat{\theta}_i) = 0

(11)

となり、各ユニット $i$ について

\sum^T_{t=1} (Y_{it} - \hat{\beta}_0 - \hat{\beta}_1 X_{it} - \hat{\theta}_i) = 0\\ \to \sum^T_{t=1} (Y_{it} - \hat{\beta}_0 - \hat{\beta}_1 X_{it}) = T \hat{\theta}_i\\ \begin{align} \to \hat{\theta}_i &= \frac{1}{T} \sum^T_{t=1} (Y_{it} - \hat{\beta}_0 - \hat{\beta}_1 X_{it})\\ &= \frac{1}{T} \sum^T_{t=1} Y_{it} - \hat{\beta}_0 - \hat{\beta}_1 \frac{1}{T} \sum^T_{t=1} X_{it}\\ &= \bar{Y}_{i} - \hat{\beta}_0 - \hat{\beta}_1 \bar{X}_{i}, \hspace{1em} \left(\bar{Y}_{i} := \frac{1}{T} \sum^T_{t=1} Y_{it}, \bar{X}_{i} := \frac{1}{T} \sum^T_{t=1} X_{it}\right)\\ \end{align}

(12)

となる。これを一階条件に代入すると

\begin{align} \sum^N_{i=1} \sum^T_{t=1} X_{it} (Y_{it} - \hat{\beta}_0 - \hat{\beta}_1 X_{it} - \bar{Y}_{i} + \hat{\beta}_0 + \hat{\beta}_1 \bar{X}_{i}) &= 0 \\ \sum^N_{i=1} \sum^T_{t=1} X_{it} (Y_{it} - \bar{Y}_{i} - \hat{\beta}_1 X_{it} + \hat{\beta}_1 \bar{X}_{i}) &= 0 \\ \sum^N_{i=1} \sum^T_{t=1} X_{it} (Y_{it} - \bar{Y}_{i} - \hat{\beta}_1 (X_{it} - \bar{X}_{i})) &= 0 \\ \sum^N_{i=1} \sum^T_{t=1} X_{it} (Y_{it} - \bar{Y}_{i}) - \sum^N_{i=1} \sum^T_{t=1} X_{it} \hat{\beta}_1 (X_{it} - \bar{X}_{i}) &= 0 \end{align}

(13)

\sum^N_{i=1} \sum^T_{t=1} X_{it} \hat{\beta}_1 (X_{it} - \bar{X}_{i}) = \sum^N_{i=1} \sum^T_{t=1} X_{it} (Y_{it} - \bar{Y}_{i})

(14)

こう移行できるか？

\frac{1}{NT} \sum^N_{i=1} \sum^T_{t=1} X_{it} \hat{\beta}_1 (X_{it} - \bar{X}_{i}) = \frac{1}{NT} \sum^N_{i=1} \sum^T_{t=1} X_{it} (Y_{it} - \bar{Y}_{i}) \\ \hat{\beta}_1 \frac{1}{NT} \sum^N_{i=1} \sum^T_{t=1} X_{it} (X_{it} - \bar{X}_{i}) = \frac{1}{NT} \sum^N_{i=1} \sum^T_{t=1} X_{it} (Y_{it} - \bar{Y}_{i})

(15)

\hat{\beta}_1 = \frac{ \sum^N_{i=1} \sum^T_{t=1} X_{it} (Y_{it} - \bar{Y}_{i}) }{ \sum^N_{i=1} \sum^T_{t=1} X_{it} (X_{it} - \bar{X}_{i}) }

(16)

TODO: このへん整理

y = data.value
x = data.inv
sum(x * (y - y.mean())) / sum(x * (x - x.mean()))

5.193824955195333

two-way effect model¶

時点効果＋個体効果といった2つの効果を同時に固定する

Y_{it} =\beta_1 X_{it} +\theta_i + \pi_t + \varepsilon_{it}

(17)

個体の固定効果 $\theta_i$ と時間の固定効果 $\pi_t$ の両方を除去したい場合は、それぞれの推定方法の組み合わせになる。

$N-1$ 個の個体ダミー変数と $T-1$ 個の時間ダミー変数を用いたOLS推定
entity demeaningと $T-1$ 個の時間ダミー変数を用いたOLS推定
time demeaningと $N-1$ 個の個体ダミー変数を用いたOLS推定
entity & time demeaningを用いたOLS推定
- 説明変数と被説明変数について、個体と時間両方の平均を引いてOLS推定

なお、パネルデータを活用した計量経済分析では、時間固定効果がないと仮定できるケースはまれであるため、通常はone-way固定効果モデルではなくtwo-way固定効果モデルを用いる。

data = sm.datasets.get_rdataset("Grunfeld", package="plm").data

rows = []
for group in ["year", "firm"]:
    for _, d in data.groupby(group):
        for col in ["value", "inv", "capital"]:
            d[col] = (d[col] - d[col].mean())
    rows.append(d)
df = pd.concat(rows)

two_way_entity_demeaned = sm.OLS.from_formula(formula="inv ~ -1 + value + capital", data=df).fit()
two_way_entity_demeaned.summary()

Loading...

比較¶

from stargazer.stargazer import Stargazer

sg = Stargazer([pooled_ols, first_diff, lsdv, entity_demeaned])
sg.custom_columns(labels=["pooled", "first_diff", "lsdv", "entity_demeaned"], separators=[1] * 4)
sg

Loading...

変量効果モデル¶

変量効果という概念もある

個体に固有の効果が説明変数と相関する場合、その効果を**固定効果（fixed effect）**と呼ぶ。
個体に固有の効果が説明変数と相関しない場合、その効果を**変量効果（random effect）**と呼ぶ。

Pooled OLSでも推定できる¶

$\theta_i$ を個体に固有の効果として

Y_{it} = \boldsymbol{\beta} X_{it} + \theta_i + \varepsilon_{it}

(18)

という線形モデルを考えたとき、変量効果モデルでは $\theta_i$ は $X_{it}$ と無相関であるため $\theta_i + \varepsilon_{it}$ をひとまとめに誤差項 $\epsilon_{it}$ として捉えて

Y_{it} = \boldsymbol{\beta} X_{it} + \epsilon_{it}

(19)

のように扱いPooled OLSとして推定すると、通常の最小二乗法の仮定を満たすためPooled OLS推定量は一致性をもつ。

変量効果モデル¶

通常、変量効果モデルと呼ぶ場合は $\theta_i$ と $X_{it}$ の独立を仮定し、FGLSで推定を行うらしい

参考文献¶

mixtape
西山ほか（2019）『計量経済学』有斐閣
藪友良パネル分析

パネルデータ¶

Pooled OLS¶

固定効果モデル¶

one-way fixed effect model¶

個体固定効果モデル¶

(1) "一回の階差モデル（first difference model）"によるOLS推定¶

(2) “N−1N-1N−1個のダミー説明変数”を用いたOLS推定¶

(3) ”平均差分法（Entity-demeaned）”を用いたOLS推定¶

導出¶

two-way effect model¶

比較¶

変量効果モデル¶

Pooled OLSでも推定できる¶

変量効果モデル¶

参考文献¶

(2) “ $N-1$ 個のダミー説明変数”を用いたOLS推定¶