固定効果モデル#

パネルデータ#

異なる時点で同じ個体についての情報を観測したデータをパネルデータ (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 \]

のような形で表記される

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: 資本ストック
firm year inv value capital
0 1 1935 317.60 3078.50 2.80
1 1 1936 391.80 4661.70 52.60
2 1 1937 410.60 5387.10 156.90
3 1 1938 257.70 2792.20 209.20
4 1 1939 330.80 4313.20 203.40
... ... ... ... ... ...
195 10 1950 3.42 69.05 8.74
196 10 1951 4.67 83.04 9.07
197 10 1952 6.00 74.42 9.93
198 10 1953 6.53 63.51 11.68
199 10 1954 5.12 58.12 14.33

200 rows × 5 columns

Pooled OLS#

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

\[ Y_{it} = \boldsymbol{\beta} X_{it} + \varepsilon_{it} \]
pooled_ols = sm.OLS.from_formula(formula="inv ~ value + capital", data=data).fit()
pooled_ols.summary()
OLS Regression Results
Dep. Variable: inv R-squared: 0.812
Model: OLS Adj. R-squared: 0.811
Method: Least Squares F-statistic: 426.6
Date: Wed, 08 May 2024 Prob (F-statistic): 2.58e-72
Time: 23:28:52 Log-Likelihood: -1191.8
No. Observations: 200 AIC: 2390.
Df Residuals: 197 BIC: 2399.
Df Model: 2
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept -42.7144 9.512 -4.491 0.000 -61.472 -23.957
value 0.1156 0.006 19.803 0.000 0.104 0.127
capital 0.2307 0.025 9.055 0.000 0.180 0.281
Omnibus: 27.240 Durbin-Watson: 0.358
Prob(Omnibus): 0.000 Jarque-Bera (JB): 89.572
Skew: 0.469 Prob(JB): 3.55e-20
Kurtosis: 6.141 Cond. No. 2.46e+03


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 2.46e+03. This might indicate that there are
strong multicollinearity or other numerical problems.

固定効果モデル#

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

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

個体内の変動を使うため、固定効果推定量は「群内(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) \]

パネルデータを用いることができる場合、以下の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}) \]

記号を置き換えて、

\[ \Delta Y_{it} = \beta_1 \Delta X_{it}+ \Delta \varepsilon_{it} \]
  • 推定方法:

    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
firm year inv value capital
0 1 1936 74.20 1583.20 49.80
1 1 1954 182.30 -648.10 449.00
2 1 1953 413.20 1316.80 346.80
3 1 1952 135.30 91.90 222.80
4 1 1951 113.00 1077.40 108.70
... ... ... ... ... ...
185 10 1937 0.19 -5.74 -0.14
186 10 1936 -0.54 17.03 0.21
187 10 1953 0.53 -10.91 1.75
188 10 1944 0.25 -0.42 -0.17
189 10 1954 -1.41 -5.39 2.65

190 rows × 5 columns

first_diff = sm.OLS.from_formula(formula="inv ~ -1 + value + capital", data=delta).fit()
first_diff.summary()
OLS Regression Results
Dep. Variable: inv R-squared (uncentered): 0.429
Model: OLS Adj. R-squared (uncentered): 0.423
Method: Least Squares F-statistic: 70.58
Date: Wed, 08 May 2024 Prob (F-statistic): 1.36e-23
Time: 23:28:52 Log-Likelihood: -982.76
No. Observations: 190 AIC: 1970.
Df Residuals: 188 BIC: 1976.
Df Model: 2
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
value 0.0891 0.008 10.816 0.000 0.073 0.105
capital 0.2787 0.047 5.910 0.000 0.186 0.372
Omnibus: 41.981 Durbin-Watson: 1.782
Prob(Omnibus): 0.000 Jarque-Bera (JB): 403.556
Skew: 0.398 Prob(JB): 2.34e-88
Kurtosis: 10.095 Cond. No. 5.75


Notes:
[1] R² is computed without centering (uncentered) since the model does not contain a constant.
[2] Standard Errors assume that the covariance matrix of the errors is correctly specified.

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

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

\[\begin{split} 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.} \end{split}\]
  • 推定方法:

    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()
OLS Regression Results
Dep. Variable: inv R-squared: 0.944
Model: OLS Adj. R-squared: 0.941
Method: Least Squares F-statistic: 288.5
Date: Wed, 08 May 2024 Prob (F-statistic): 2.41e-111
Time: 23:28:52 Log-Likelihood: -1070.8
No. Observations: 200 AIC: 2166.
Df Residuals: 188 BIC: 2205.
Df Model: 11
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept -70.2967 49.708 -1.414 0.159 -168.354 27.760
firm[T.2] 172.2025 31.161 5.526 0.000 110.732 233.673
firm[T.3] -165.2751 31.776 -5.201 0.000 -227.958 -102.593
firm[T.4] 42.4874 43.910 0.968 0.334 -44.132 129.107
firm[T.5] -44.3201 50.492 -0.878 0.381 -143.924 55.284
firm[T.6] 47.1354 46.811 1.007 0.315 -45.206 139.477
firm[T.7] 3.7432 50.565 0.074 0.941 -96.004 103.491
firm[T.8] 12.7511 44.053 0.289 0.773 -74.150 99.652
firm[T.9] -16.9256 48.453 -0.349 0.727 -112.508 78.656
firm[T.10] 63.7289 50.330 1.266 0.207 -35.556 163.013
value 0.1101 0.012 9.288 0.000 0.087 0.134
capital 0.3101 0.017 17.867 0.000 0.276 0.344
Omnibus: 29.604 Durbin-Watson: 1.079
Prob(Omnibus): 0.000 Jarque-Bera (JB): 162.947
Skew: 0.283 Prob(JB): 4.13e-36
Kurtosis: 7.386 Cond. No. 6.42e+04


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 6.42e+04. This might indicate that there are
strong multicollinearity or other numerical problems.

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

\[\begin{split} \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} \end{split}\]
  • 推定方法:

    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()
OLS Regression Results
Dep. Variable: inv R-squared (uncentered): 0.767
Model: OLS Adj. R-squared (uncentered): 0.764
Method: Least Squares F-statistic: 325.5
Date: Wed, 08 May 2024 Prob (F-statistic): 2.59e-63
Time: 23:28:53 Log-Likelihood: -1070.8
No. Observations: 200 AIC: 2146.
Df Residuals: 198 BIC: 2152.
Df Model: 2
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
value 0.1101 0.012 9.532 0.000 0.087 0.133
capital 0.3101 0.017 18.336 0.000 0.277 0.343
Omnibus: 29.604 Durbin-Watson: 1.079
Prob(Omnibus): 0.000 Jarque-Bera (JB): 162.947
Skew: 0.283 Prob(JB): 4.13e-36
Kurtosis: 7.386 Cond. No. 1.74


Notes:
[1] R² is computed without centering (uncentered) since the model does not contain a constant.
[2] Standard Errors assume that the covariance matrix of the errors is correctly specified.
md = smf.mixedlm("inv ~ value + capital", data, groups=data["year"])
mdf = md.fit()
print(mdf.summary())
/usr/local/lib/python3.10/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.10/site-packages/statsmodels/regression/mixed_linear_model.py:2200: ConvergenceWarning: Retrying MixedLM optimization with lbfgs
  warnings.warn(
/usr/local/lib/python3.10/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.10/site-packages/statsmodels/regression/mixed_linear_model.py:2200: 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.10/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.10/site-packages/statsmodels/regression/mixed_linear_model.py:2206: ConvergenceWarning: MixedLM optimization failed, trying a different optimizer may help.
  warnings.warn(msg, ConvergenceWarning)
/usr/local/lib/python3.10/site-packages/statsmodels/regression/mixed_linear_model.py:2218: ConvergenceWarning: Gradient optimization failed, |grad| = 1.131047
  warnings.warn(msg, ConvergenceWarning)
/usr/local/lib/python3.10/site-packages/statsmodels/regression/mixed_linear_model.py:2261: 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} \]

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

\[ (\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 \]

となり、各ユニット\(i\)について

\[\begin{split} \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} \end{split}\]

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

\[\begin{split} \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} \end{split}\]
\[ \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}) \]

こう移行できるか?

\[\begin{split} \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}) \end{split}\]
\[ \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}) } \]

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} \]

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

  1. \(N-1\)個の個体ダミー変数と\(T-1\)個の時間ダミー変数を用いたOLS推定

  2. entity demeaningと\(T-1\)個の時間ダミー変数を用いたOLS推定

  3. time demeaningと\(N-1\)個の個体ダミー変数を用いたOLS推定

  4. 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()
OLS Regression Results
Dep. Variable: inv R-squared (uncentered): 0.898
Model: OLS Adj. R-squared (uncentered): 0.891
Method: Least Squares F-statistic: 123.8
Date: Wed, 08 May 2024 Prob (F-statistic): 1.24e-14
Time: 23:28:53 Log-Likelihood: -173.11
No. Observations: 30 AIC: 350.2
Df Residuals: 28 BIC: 353.0
Df Model: 2
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
value 0.1345 0.035 3.818 0.001 0.062 0.207
capital 0.3315 0.097 3.431 0.002 0.134 0.529
Omnibus: 39.780 Durbin-Watson: 2.574
Prob(Omnibus): 0.000 Jarque-Bera (JB): 167.640
Skew: -2.559 Prob(JB): 3.96e-37
Kurtosis: 13.389 Cond. No. 6.79


Notes:
[1] R² is computed without centering (uncentered) since the model does not contain a constant.
[2] Standard Errors assume that the covariance matrix of the errors is correctly specified.

比較#

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
Dependent variable: inv
pooledfirst_difflsdventity_demeaned
(1)(2)(3)(4)
Intercept-42.714***-70.297
(9.512)(49.708)
capital0.231***0.279***0.310***0.310***
(0.025)(0.047)(0.017)(0.017)
firm[T.10]63.729
(50.330)
firm[T.2]172.203***
(31.161)
firm[T.3]-165.275***
(31.776)
firm[T.4]42.487
(43.910)
firm[T.5]-44.320
(50.492)
firm[T.6]47.135
(46.811)
firm[T.7]3.743
(50.565)
firm[T.8]12.751
(44.053)
firm[T.9]-16.926
(48.453)
value0.116***0.089***0.110***0.110***
(0.006)(0.008)(0.012)(0.012)
Observations200190200200
R20.8120.4290.9440.767
Adjusted R20.8110.4230.9410.764
Residual Std. Error94.408 (df=197)42.896 (df=188)52.768 (df=188)51.418 (df=198)
F Statistic426.576*** (df=2; 197)70.578*** (df=2; 188)288.500*** (df=11; 188)325.451*** (df=2; 198)
Note:*p<0.1; **p<0.05; ***p<0.01

変量効果モデル#

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

  • 個体に固有の効果が説明変数と相関する場合、その効果を**固定効果(fixed effect)**と呼ぶ。

  • 個体に固有の効果が説明変数と相関しない場合、その効果を**変量効果(random effect)**と呼ぶ。

Pooled OLSでも推定できる#

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

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

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

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

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

変量効果モデル#

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

参考文献#