17_Panel.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# パネル・データ分析"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div name=\"html-admonition\" style=\"font-size: 0.8em\">\n",
    "<input type=\"button\" onclick=\"location.href='https://translate.google.com/translate?hl=&sl=ja&tl=en&u='+window.location;\" value=\"Google translation\" style=\"color:#ffffff;background-color:#008080; height:25px\" onmouseover=\"this.style.background='#99ccff'\" onmouseout=\"this.style.background='#008080'\"/> in English or the language of your choice.\n",
    "</div><br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.formula.api as smf\n",
    "import wooldridge\n",
    "\n",
    "from collections import OrderedDict\n",
    "from linearmodels.panel.data import PanelData\n",
    "from linearmodels.panel import PanelOLS, PooledOLS, RandomEffects, compare\n",
    "from scipy.stats import gaussian_kde, norm, t\n",
    "\n",
    "# 警告メッセージを非表示\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "パネル・データを使った次のモデルについて説明する。\n",
    "* 固定効果モデル（Fixed Effects Model）\n",
    "* ランダム効果モデル（Random Effects Model）\n",
    "* 相関ランダム効果モデル（Correlated Random Effects Model） "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "## 固定効果モデル"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true,
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "### 説明"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "パネル・データを使う場合の問題は観察単位の異質性を捉える変数 $a_i$ が説明変数 $x_{it}$ と相関することにより，単純なOLS推定量に異質性バイアスが発生することである。その対処方法として，回帰式から$a_i$をなくす１階差分推定について説明した。この章では，同じように$a_i$をなくす代替推定法として固定効果推定（Fixed Effects Estimator）について考える。まず，その基礎となる固定効果モデル（Fixed Effects Model; FEモデル）について解説する。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "次式を考えよう。\n",
    "\n",
    "$$y_{it}= \\beta_0 + a_i + \\beta_1x_{it}+u_{it}\\qquad i=1,2,..,n\\quad t=0,1,...,T\\qquad\\qquad\\left(\\text{式１}\\right)$$\n",
    "\n",
    "$a_i$は時間に依存せず，個体$i$特有の因子となっている。この$i$に固定された$a_i$が重要なモデルの要素なっているため固定効果モデルと呼ばれる。\n",
    "\n",
    "$i$を所与として，両辺のそれぞれの変数の時間に対しての平均を計算すると次式となる。\n",
    "\n",
    "$$\\bar{y}_i= \\beta_0 + a_i + \\beta_1\\bar{x}_i+\\bar{u}_i\\qquad\\qquad\\left(\\text{式２}\\right)$$\n",
    "\n",
    "（注意）$a_i$は時間に対して不変なため，$a_i$はそのままの値を取る。\n",
    "\n",
    "(式１)と(式２)の差を取ると次の固定効果推定式が導出できる。\n",
    "\n",
    "$$\\ddot{y}_i = \\beta_1\\ddot{x}_i+\\ddot{u}_i\\qquad\\qquad\\left(\\text{式３}\\right)$$\n",
    "\n",
    "ここで\n",
    "* $\\ddot{y}_i=y_i-\\bar{y}$\n",
    "* $\\ddot{x}_i=x_i-\\bar{x}$\n",
    "* $\\ddot{u}_i=u_i-\\bar{u}$\n",
    "\n",
    "は平均からの乖離（demeaned values）。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "＜仮定と推定量の性質＞\n",
    "1. (式１)のように回帰式は線形\n",
    "1. 標本のランダム抽出\n",
    "1. $x_{it}$は時間と共に変動し、説明変数間の完全多重共線性は存在しない。\n",
    "1. $\\text{E}(u_{it}|X_i,a_i)=0$、即ち、$a_i$ と $x_{it}$, $t=1,2,..,T$ を所与として $u_{it}$ の平均は0\n",
    "\n",
    "この仮定の下で，OLS推定量 $\\hat{\\beta}_1$ は\n",
    "* 不偏性を満たす。\n",
    "* $T$を一定として標本の大きさが十分に大きい場合，一致性を満たす。\n",
    "* 固定効果推定量（Fixed Effects Estimator or FE Estimator）もしくは Within Estimatorと呼ばれる。\n",
    "\n",
    "\n",
    "<br>\n",
    "\n",
    "---\n",
    "＜良い点＞\n",
    "* 推定量は上述の仮定の下で不偏性・一致性が満たされる。\n",
    "\n",
    "＜悪い点＞\n",
    "* 時間に対して一定の説明変数は，$a_i$と同じように，（式３）に残らない。従って，時間不変の変数の効果を推定することはできない。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 「手計算」による推定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`wagepan`というデータセットを使い，賃金に対する労働組合などの変数の効果を推定する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true,
    "tags": [
     "output_scroll"
    ]
   },
   "outputs": [],
   "source": [
    "wagepan = wooldridge.data('wagepan')\n",
    "wooldridge.data('wagepan', description=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "被説明変数\n",
    "* `lwage`：賃金（対数）\n",
    "\n",
    "説明変数\n",
    "* `union`：労働組合参加（ダミー変数）\n",
    "* `married`：未・既婚（ダミー変数）\n",
    "* `exper`：労働市場参加年数\n",
    "*  `d81`：1981年のダミー変数\n",
    "*  `d82`：1982年のダミー変数\n",
    "*  `d83`：1983年のダミー変数\n",
    "*  `d84`：1984年のダミー変数\n",
    "*  `d85`：1985年のダミー変数\n",
    "*  `d86`：1986年のダミー変数\n",
    "*  `d87`：1987年のダミー変数\n",
    "\n",
    "`DataFrame`の`groupby`を使って変数をグループ化する際に使う変数\n",
    "* `nr`：労働者のID"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "```{admonition} コメント\n",
    ":class: tip\n",
    "時間に対して変化しない変数は使えない。例えば，\n",
    "* `educ`（説明変数に入れば教育の収益率が推定可能である。）\n",
    "* `black`，`hisp`（を使うと人種間の賃金格差も推定できる。）\n",
    "\n",
    "は固定効果モデル回帰式には入れることができない。\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# 説明変数のリスト\n",
    "exog = ['married','union','expersq','d81','d82','d83','d84','d85','d86','d87']\n",
    "\n",
    "# 全ての変数のリスト\n",
    "var = ['lwage']+exog\n",
    "\n",
    "# 使う変数だけで構成されるDataFrame\n",
    "df = wagepan.loc[:,['nr']+var]\n",
    "\n",
    "# varの平均からの乖離を計算（下の説明（１）を参照）\n",
    "df_g = df.groupby('nr')\n",
    "df_mean = df_g[var].transform('mean')\n",
    "df_md = df.loc[:,var]-df_mean\n",
    "\n",
    "# 説明変数の行列（下の説明（２）を参照）\n",
    "X = df_md.loc[:,exog].values\n",
    "\n",
    "# 被説明変数のベクトル\n",
    "Y = df_md.loc[:,'lwage'].values\n",
    "\n",
    "# OLSの計算\n",
    "params = np.linalg.inv((X.T)@X)@(X.T)@Y\n",
    "\n",
    "# 結果の表示（下の説明（３）を参照）\n",
    "for idx, name in enumerate(exog):\n",
    "    print(f'{name}: {params[idx]:.4}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "説明（１）\n",
    "* `df.groupby('nr')`：`nr`でグループ化したオブジェクトを作成する（番外編のGapminderを参照）。\n",
    "* `df_g[var].transform('mean')`\n",
    "    * `df_g[var]`：グループ化計算のために`var`の変数だけを使うことを指定する。\n",
    "    * `.transform('mean')`：指定された変数（`var`）のグループ内平均（それぞれの`nr`内での平均）を計算し，その平均で構成される`DataFrame`を作成する。作成された`DataFrame`の行数は`df`と同じになり，グループ内平均が同じ列内でリピートされることになる。\n",
    "* `df.loc[:,var]-df_mean`：`var`のそれぞれの変数の平均からの乖離を計算する。\n",
    "* `df_mean = ...`と`df_md = ...`の2行を次のように1行にまとめることも可能。\n",
    "    * `df_md = df_g[var].transform(lambda x: x-x.mean())`\n",
    "\n",
    "説明（２）\n",
    "* `.values`は`DataFrame`を`numpy`の`array`として返す。\n",
    "\n",
    "説明（３）\n",
    "* `enumerate(exog)`は`exog`の変数のインデックスと変数名の両方をタプルとして返す。\n",
    "    * `for`ループで`enumerate()`は，引数にリストや`array`を入れると2つの値を返す。\n",
    "    * 1つ目の返り値：要素のインデックス\n",
    "    * 2つ目の返り値：要素自体\n",
    "    * 例えば，`enumerate(['A','B','C'])`であれば，返り値は`(0,'A')`，`(1,'B')`，`(2,'C')`となる。\n",
    "* `f'{name}: {params[idx]:.4}'`\n",
    "    * `f`：以前説明したf-string\n",
    "    * `{name}`：`exog`の変数名を代入する。\n",
    "    * `{params[idx]:.4}`\n",
    "        * `params[idx]`：`idx`番目の`params`を代入する。\n",
    "        * `:.4`：小数点第4位までを表示する。\n",
    "\n",
    "---\n",
    "「手計算」の結果は`linearmodels`を使う結果と同じになることを下で確認する。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### `linearmodels`を使う計算"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### EntityEffects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`linearmodels`には`PanelOLS`モジュールがあり，その関数`from_formula()`を使うことにより，`statsmodels`同様，回帰式を`y ~ x`の文字列で書くことが可能となる。その際，次の点に注意すること。\n",
    "* 固定効果推定をする場合，回帰式に`+ EntityEffect`を含める。\n",
    "    * このオプションにより変数の平均からの乖離は自動で計算されることになる。\n",
    "    * `+ EntityEffect`を含めないと`PooledOLS`（通常のOLS）と等しくなる。\n",
    "\n",
    "\n",
    "まず`wagepan`を`MultiIndex`化する。これにより`linearmodels`を使いFD推定可能となる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "wagepan = wagepan.set_index(['nr','year'],drop=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true,
    "tags": [
     "output_scroll"
    ]
   },
   "outputs": [],
   "source": [
    "wagepan.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true,
    "tags": [
     "output_scroll"
    ]
   },
   "outputs": [],
   "source": [
    "wagepan.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "次に`PanelData`オブジェクトに変換しデータの特徴を調べる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "wagepanp = PanelData(wagepan)\n",
    "wagepanp.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "* 44: 変数の数\n",
    "* 8: 期間数（年）\n",
    "* 545：観察単位の数（人数）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "次に，balanced もしくは unbalanced data set かを確認する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "(~wagepanp.isnull).all()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "このデータ・セットはbalancedだが，unbalanced だったとしても，固定効果モデルの考え方や以下で説明するコードは変わらない。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "実際に回帰式を書くことにする。使い方は`statsmodels`と似ている。\n",
    "* `PanelOLS`モジュールの関数`.from_formula`を使い次のように引数を指定する。\n",
    "\n",
    "$$\\text{.from_formula}(\\text{回帰式}, \\text{データ})$$\n",
    "\n",
    "* `EntityEffects`を加える。\n",
    "* 定数項を入れたい場合は，`1`を回帰式に追加する。入れなければ定数項なしの推定となる。\n",
    "\n",
    "* 以下では時間ダミー`C(year)`が入るので入れない。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "formula_fe = 'lwage ~ married + union + expersq \\\n",
    "                      +d81+d82+d83+d84+d85+d86+d87 + EntityEffects'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "* 固定効果モデルのインスタンスの作成"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "mod_fe = PanelOLS.from_formula(formula_fe, data=wagepan)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "* `statsmodels`と同じように，そこから得た結果にメソッド`.fit()`を使い計算し結果が返される。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "result_fe = mod_fe.fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "＜結果の表示方法＞\n",
    "1. `res_fe`を実行。\n",
    "1. `res_fe`に関数`print()`を使うと見やすい。\n",
    "1. `res_fe`には属性`summary`が用意されているが，表示方法1と同じ。\n",
    "1. `summary`には属性`tables`があり，２つの表がリストとして格納されている。\n",
    "    * `tables[0]`：検定統計量の表（`print()`を使うと見やすくなる）\n",
    "    * `tables[1]`：係数の推定値やp値などの表（`print()`を使うと見やすくなる）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "print(result_fe.summary.tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "（結果）\n",
    "* `exper**2`の係数が負で統計的有意性が高いのは，賃金に対して経験の効果は低減することを示している。\n",
    "* `married`の係数は正であり，優位性が全くないわけではない。賃金の既婚プレミアムと呼ばれるものである。\n",
    "* `union`は労働組合の影響を示しているが，予測通りである。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "$R^2$を表示してみる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "print(result_fe.summary.tables[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "＜この表にある$R^2$について＞\n",
    "* $R^2$＝$R^2(\\text{Within})$：(式３)を推定した際の$R^2$である（$\\ddot{y}_i$が$\\ddot{x}_i$によってどれだけ説明されたかを示す）。\n",
    "* $R^2(\\text{Between})$：(式２)を推定した際の$R^2$である（$\\hat{y}_i$が$\\hat{x}_i$によってどれだけ説明されたかを示す）。\n",
    "* $R^2(\\text{Overall})$：(式１)を推定した際の$R^2$である（${y}_i$が${x}_i$によってどれだけ説明されたかを示す）。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "````{hint}\n",
    "(式３)の回帰式には定数項はないため，上の計算で使った回帰式`formula_fe`には定数項を表す`1`は省かれている。(式３)と整合性が取れたコードとなっており，直観的にも分かりやすい。では，`formula_fe`に`1`を加えて次のように書き換えて計算するとどうなるのだろうか。\n",
    "```\n",
    "formula_fe = 'lwage ~ 1 + married + union + expersq \\\n",
    "                      +d81+d82+d83+d84+d85+d86+d87 + EntityEffects'\n",
    "```\n",
    "このコードを使うと定数項の推定値が計算されることになるが，(式３)とは整合性がないように見える。実際そうなのだが，`linearmodels`では，この場合の定数項の推定値は，(式３)の個体固有の因子である$a_i$の平均値を表す仕様になっている。`formula_fe`に`1+`が追加されても，他の変数の推定値などは変わらない。確かめてみよう！\n",
    "````"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### TimeEffects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "上の推定式では時間ダミー変数として使い，観察単位全てに共通な時間的な影響を捉えた。具体的には，インフレにより賃金は変化するが，その変化は全ての労働者には対して同じであり，その効果を時間ダミー変数が捉えている。それを**時間効果**と呼ぶ。このような時間ダミー変数を加えた理由は，時間効果が他の変数（例えば，`married`）の係数を「汚さない」ようにするためであり，よりピュアの効果を推定するためである。一方，`linearmodels`では，わざわざ時間ダミー変数を作らずとも`TimeEffects`を回帰式に追加することにより，時間効果を自動的に「排除」することができる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "formula_fe2 = 'lwage ~ married + union + expersq + TimeEffects + EntityEffects'\n",
    "result_fe2 = PanelOLS.from_formula(formula_fe2, data=wagepan).fit()\n",
    "print(result_fe2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`result_fe`と同じ結果を確認できる。$R^2$の値は少し変わっているが，これは時間ダミーを入れて計算している訳ではないためである。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "## ダミー変数モデルとしての固定効果モデル"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "固定効果推定量は，他の推定方法でも計算することができる。その１つがダミー変数推定である。アイデアは簡単で，観察単位のダミー変数を使い異質性を捉えるのである。推定方法も簡単で，（式１）に観察単位のダミー変数を加えて通常のOLS推定をおこなうだけである。\n",
    "\n",
    "上で使った`wagepan`を使い推定する。まず，観察単位のダミー変数として推定式に`C(nr)`を追加する"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "formula_dum = 'lwage ~  1 + married + union + expersq \\\n",
    "                        +d81+d82+d83+d84+d85+d86+d87 + C(nr)'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`PooledOLS`モジュールの関数`from_formula`を使って推定式を定義する。ここで`PooledOLS`とは，`statsmodels`で使う通常のOLS推定と同じである。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "result_dum = PooledOLS.from_formula(formula_dum, data=wagepan).fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`nr`のダミー変数が544あるため，そのまま結果を表示せずに，以下のセルにあるフィルターを使いメインの変数だけを表示する。下のコードの`filter`で使う（`result_dum`に続く）属性・メソッドについて：\n",
    "* `params`：パラメーターを取得する属性\n",
    "* `index`：`params`のインデックスを取得する属性\n",
    "* `str`：`index`にある文字列を操作可能なオブジェクトとして取得する属性\n",
    "* `contains()`：引数の文字列が文字列の中にある場合`True`を返すメソッド\n",
    "* `( )`：括弧の中を先に評価するという意味\n",
    "* `tolist()`：リストに変換するメソッド"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "filter = (result_dum.params.index.str.contains('nr') == False).tolist()\n",
    "result_dum.params[filter]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "３つの変数`I(exper**2)`，`married`，`union`の係数は固定効果モデルと等しいことが確認できる。\n",
    "\n",
    "t値とp値の表示には次のコードを使う。\n",
    "* t値：`result_dum.tstats[filter]`\n",
    "* p値：`result_dum.pvalues[filter]`\n",
    "\n",
    "この場合，`educ`等の時間に対して不変の変数も推定式に加えることが可能である。一方，この方法はパラメーターの数は544+11=555あり，自由度が非常に低くなるのが欠点である。\n",
    "\n",
    "\n",
    "（注意）\n",
    "\n",
    "この推定式に時間に対して不変の変数（例えば，`educ`，`blac`，`hisp`，`exper`）を追加すると推定できない（エラーが発生する）。理由は，それらの変数は`nr`ダミー変数と完全に同じ動きをするためである。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "## ランダム効果モデル"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true,
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "### 説明"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "ここでは，パネル・データを使い推定するランダム効果モデル（Random Effects Model; RE Model）を解説する。もう一度（式１）を考えよう。\n",
    "\n",
    "（仮定） $a_i$は観察単位 $i$ に対して一定であるが，$i$によってその値は異なる。\n",
    "\n",
    "＜固定効果モデルの仮定＞\n",
    "* $\\text{Cov}\\left(a_ix_{it}\\right)\\neq 0,\\quad t=1,2,...,T$\n",
    "\n",
    "＜ランダム効果モデルの仮定＞\n",
    "* $\\text{Cov}\\left(a_ix_{it}\\right)=0,\\quad t=1,2,...,T$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "この違いを念頭に（式１）を次式に書き換える。\n",
    "\n",
    "$$y_{it}= \\beta_0 + \\beta_1x_{it}+e_{it}\\qquad\\qquad\\left(\\text{式４}\\right)$$\n",
    "\n",
    "ここで，\n",
    "\n",
    "$$e_{it}=a_i+u_{it}$$\n",
    "\n",
    "しかし，$e_{it}$と$e_{its},\\;t\\ne s$の共分散は\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\text{Cov}\\left(e_{it},e_{is}\\right)\n",
    "&=\\text{Cov}\\left(a_i+u_{it},a_i+u_{is}\\right)\\\\\n",
    "&=\\text{Cov}\\left(a_i,a_i\\right)+\\text{Cov}\\left(a_i,u_{is}\\right)\n",
    "+\\text{Cov}\\left(u_{it},a_i\\right)+\\text{Cov}\\left(u_{it},u_{is}\\right)\\\\\n",
    "&=\\text{Cov}\\left(a_{i},a_{i}\\right)+0+0+0\\\\\n",
    "&=\\text{Var}\\left(a_{i}\\right)\\\\\n",
    "&\\neq 0\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "となる。即ち，誤差項が自己相関することになる。自己相関係数を計算してみよう。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\rho_{e_{it}e_{is}}\n",
    "=\\frac{\\text{Cov}(e_{it}e_{is})}{\\sqrt{\\text{Var}(e_{it})\\text{Var}(e_{is})}}\n",
    "=\\frac{\\sigma_{a}^2}{\\sigma_{u}^2+\\sigma_a^2}\n",
    "$$\n",
    "\n",
    "* $\\sigma_{a}^2=\\text{Var}(a_i)$（$a_{i}$の分散）\n",
    "* $\\sigma_{u}^2=\\text{Var}(u_{it})$（$u_{it}$の分散）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "＜含意＞\n",
    "* $a_i$と$x_{it}$に相関がないため次の推定方法で一致性を満たす推定量を計算できる。\n",
    "    * １つの時点でお横断面データを使いOLS推定（でも他のデータはどうする？）\n",
    "    * 全てのデータを使い何の区別もなくプールするPooled OLS推定（しかしパネル・データの特性を有効利用していない）。\n",
    "* 時系列データを扱う場合，自己相関によりOLS推定量は効率性が低くなる。これは横断面データを扱う際の不均一分散が引き起こす問題と似ている。\n",
    "* １階差分推定や固定効果推定を使う必要はなく，使うと効率性が低い推定量となる。\n",
    "    * １階差分推定や固定効果推定は$\\text{Cov}\\left(a_ix_{it}\\right)\\neq 0,\\quad t=1,2,...,T$の場合に使う。\n",
    "* 従って，問題は次の点：\n",
    "    * より効率性が高い推定方法はどのようなものか。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "その方法がランダム効果推定と言われるもので，固定効果のように，変数を平均からの**部分的**な乖離に変換することにより可能となる。具体的には，次式が推定式となる。\n",
    "\n",
    "$$\\overset{\\circ}{y}_{it}=\\beta_0(1-\\theta)+\\beta_1\\overset{\\circ}{x}_{it}+\\overset{\\circ}{e}_{it}\\qquad\\qquad\\left(\\text{式５}\\right)$$\n",
    "\n",
    "ここで\n",
    "* $\\overset{\\circ}{y}_{it}=y_{it}-\\theta\\bar{y}_i$\n",
    "* $\\overset{\\circ}{x}_{it}=x_{it}-\\theta\\bar{x}_i$\n",
    "* $\\overset{\\circ}{e}_{it}=e_{it}-\\theta\\bar{e}_i$\n",
    "\n",
    "は変数の平均からの部分的な乖離であり，乖離の度合いを決める変数$\\theta$は\n",
    "\n",
    " $$\\theta = 1-\\sqrt{\\frac{\\sigma_u^2}{\\sigma_u^2+T\\sigma_a^2}}$$\n",
    "\n",
    "と定義される。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "（直観）\n",
    "* $\\sigma_a^2=0$の場合$\\theta=0$となり，$a_i$は$i$に対して一定であり観察単位の異質性はないということである。その場合は，通常のOLSで推定することがベストとなる。\n",
    "* $\\sigma_a^2$が増加すると$\\theta$は大きくなり，より大きな平均からの乖離が必要となる。極端な場合，$\\sigma_a^2$が無限大に近づくと，$\\theta=1$となり$\\overset{\\circ}{y}=\\ddot{y}$となる。即ち，固定効果モデルはランダム効果モデルの極端な場合と解釈できる。\n",
    "\n",
    "（注意）\n",
    "* $\\theta$は事前にはわからないため推定する必要がある（`linearmodels`が自動的に計算する）。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "＜仮定と推定量の性質＞\n",
    "1. (式１)のように回帰式は線形\n",
    "1. 標本のランダム抽出\n",
    "1. 説明変数間の完全多重共線性は存在しない。\n",
    "1. 説明変数を所与として$a_i$の平均は一定（従って、$\\text{Cov}\\left(a_ix_{it}\\right)= 0$）。\n",
    "\n",
    "この仮定の下で，RE推定量 $\\hat{\\beta}_1$ は\n",
    "* 一致性を満たす。しかし不偏性は満たさない。\n",
    "\n",
    "\n",
    "<br>\n",
    "\n",
    "---\n",
    "（良い点）\n",
    "* 時間に対して不変の説明変数があってもその係数を推定できる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true,
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "### 推定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`exper`，`educ`，`black`，`hisp`を加えて回帰式を定義する。\n",
    "\n",
    "（注意）\n",
    "\n",
    "時間ダミー変数の代わりに`TimeEffects`を使わないように。\n",
    "入れてもエラーは出ませんが，`TimeEffects`は無視されて，時間ダミーがない結果と同じになります。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "formula_re = 'lwage ~ 1 + married + union + expersq \\\n",
    "                        + exper + educ + black + hisp \\\n",
    "                        + d81 + d82 + d83 +d84 + d85 + d86 + d87'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`RandomEffects`のモジュールにある関数`from_formula`を使い計算する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "result_re = RandomEffects.from_formula(formula_re, data=wagepan).fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "結果の表示。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "print(result_re.summary.tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "結果の解釈は下でする。\n",
    "\n",
    "RE推定では$\\theta$が重要な役目を果たすが，その値は`result_re`の属性`theta`を使うことにより`DataFrame`の形で表示できる。\n",
    "* balanced panel dataの場合，`theta`の値は一意で決まる。\n",
    "* unbalanced panel dataの場合，`theta`は観察単位毎に計算される。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "result_re.theta.iloc[0,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "上で説明したが$\\theta$は$u_a$と$u_it$の分散である$\\sigma_a^2$と$\\sigma_u^2$に依存する。それらの値は，属性`variance_decomposition`を表示できる。以下の返り値の内容：\n",
    "* `Effects`：$\\sigma_a^2$\n",
    "* `Residual`：$\\sigma_u^2$\n",
    "* `Percept due to Effects`：$\\dfrac{\\sigma_a^2}{\\sigma_a^2+\\sigma_u^2}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "result_re.variance_decomposition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "## 相関ランダム効果モデル"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 説明"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "相関ランダム効果モデル（CREモデル）は，固定効果モデルとランダム効果モデルの中間的な位置にあり，両方を包含している。（式１）を考えよう。更に，観察不可能な固定効果$a_i$は説明変数と次の関係にあると仮定する。\n",
    "\n",
    "$$a_i = \\alpha + \\gamma\\bar{x}_{it} + r_i\\qquad\\text{(式５)}$$\n",
    "\n",
    "* $\\bar{x}_i=\\dfrac{1}{T}\\sum_{t=1}^Tx_{it}$は説明変数の平均\n",
    "* $\\gamma$は$a_i$と$x_{it}$の相関関係を捉える係数\n",
    "* $r_i$は説明変数$x_{it}$と相関しないと仮定，即ち，$\\text{Cov}\\left(r_i\\bar{x}_{it}\\right)=0$\n",
    "\n",
    "（式５）を（式１）に代入すると次式を得る。\n",
    "\n",
    "$$y_{it}=\\alpha+\\beta x_{it} + \\gamma\\bar{x}_i + v_{it}\\qquad\\text{(式６)}$$\n",
    "\n",
    "ここで\n",
    "\n",
    "$$v_{it}=r_i + u_{it}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "（含意）\n",
    "* $\\text{Cov}\\left(r_i,\\bar{x}_{it}\\right)=0\\;\\Rightarrow\\;\\text{Cov}\\left(v_i,\\bar{x}_{it}\\right)=0$\n",
    "* REモデルと同じ構造となっており，違いは$\\bar{x}_i$が追加されている。\n",
    "    * `linearmodels`の`RandomEffects`モジュールが使える。\n",
    "* 次の結果が成立する。\n",
    "\n",
    "    $$\\hat{\\beta}_{FE}=\\hat{\\beta}_{CRE}$$\n",
    "    \n",
    "    * $\\hat{\\beta}_{FE}$：固定効果推定量\n",
    "    * $\\hat{\\beta}_{CRE}$：相関ランダム効果推定量\n",
    "    * この結果は，時間に対して不変な変数（例えば，`black`）を**含めても**成立する"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 推定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "まず $\\bar{x}_i$ を計算し，それを`wagepan`に追加する。そのために次の関数を定義する。\n",
    "\n",
    "（解説）\n",
    "\n",
    "* (1)：関数の引数\n",
    "    * `dframe`：データフレーム\n",
    "    * `ori_col`：平均を計算したい列\n",
    "    * `new_col`：計算した平均を入れる列\n",
    "* (2)：`ori_col`をグループ化し，グループ名とグループ平均からなる辞書の作成\n",
    "    * `groupby(level=0)`：行の第１インデックスに従ってグループ化\n",
    "    * `mean()`：グループ平均の計算\n",
    "    * `to_dict()`：行の第１インデックスにあるグループ名を`key`，グループ平均を`value`にする辞書の作成\n",
    "* (3)：行の第１インデックスに該当するグループ平均が並ぶリストの作成\n",
    "    * `index.get_level_values(0)`：行の第１インデックスの値を取得\n",
    "    * `to_series()`：`series`に変換\n",
    "    * `map(dict)`：`dict`の内容に従って上の`series`の値をグループ平均に入れ替える\n",
    "    * `tolist()`：リストに変換\n",
    "* (4)：`dframe`にグループ平均が並ぶ新しい列が追加し，そのラベルを`new_col`とする\n",
    "* (5)：`DataFrame`を返す"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def add_col_mean(dframe, ori_col, new_col):  # (1)\n",
    "    \n",
    "    dict = dframe.groupby(level=0)[ori_col].mean().to_dict()  # (2)\n",
    "    mean = dframe.index.get_level_values(0).to_series().map(dict).tolist()  # (3)\n",
    "    dframe.loc[:,new_col] = mean  # (4)\n",
    "    \n",
    "    return dframe   # (5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "この関数を使い，`married`，`union`，`expersq`の平均を計算し`wagepan`に追加する。\n",
    "\n",
    "（コメント）`exper`は含めない。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "wagepan = add_col_mean(wagepan, 'married', 'married_mean')\n",
    "wagepan = add_col_mean(wagepan, 'union', 'union_mean')\n",
    "wagepan = add_col_mean(wagepan, 'expersq', 'expersq_mean')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "CRE推定と結果の表示"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "formula_cre = 'lwage ~ 1 + married + union + expersq \\\n",
    "                         + married_mean + union_mean + expersq_mean \\\n",
    "                         + d81 + d82 + d83 + d84 + d85 + d86 + d87'\n",
    "\n",
    "result_cre = RandomEffects.from_formula(formula_cre, data=wagepan).fit()\n",
    "\n",
    "print(result_cre)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### ２つの利点"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "CREモデルの２つの利点を解説する。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### FE対RE検定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "FEモデルとREモデルのどちらが適しているかを調べることができるHausman検定というものがある。CREモデルを使うことにより，同様の検定が簡単に行える。（式６）を考えよう。\n",
    "* $\\gamma=0$の場合，REモデルの（式５）と同じになり，REモデルが妥当ということになる。\n",
    "* $\\gamma\\neq 0$の場合，$a_i$と$\\bar{x}_{i}$は相関することになり，これは$\\text{Cov}\\left(a_i,{x}_{it}\\right)\\neq 0$を意味し，FEモデルが妥当な推定方法となる。\n",
    "\n",
    "この考えを利用して，次の帰無仮説と対立仮説のもとで$\\gamma$の優位性を調べる。\n",
    "* $\\text{H}_0:\\;\\text{Cov}\\left(a_i,x_{it}\\right)=0$\n",
    "* $\\text{H}_a:\\;\\text{Cov}\\left(a_i,x_{it}\\right)\\neq 0$\n",
    "\n",
    "（コメント）\n",
    "* $\\gamma=0$を棄却できれば上の$\\text{H}_0$を棄却できる。\n",
    "* 平均の変数が$k$ある場合は，$\\gamma_1=\\gamma_2=\\cdots=\\gamma_k=0$を検定する。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`result_cre`のメソッド`wald_test()`を使う。引数には文字列で指定する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# 検定する係数の値を設定する\n",
    "restriction = 'married_mean = union_mean = expersq_mean = 0'\n",
    "\n",
    "# 検定結果を表示する\n",
    "result_cre.wald_test(formula=restriction)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "p値が非常に小さいので，帰無仮説は棄却できる。従って，FEモデルが妥当だと結論づけることができる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### 一定な変数を含める"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "（式６）に時間に対して一定は変数（$z_{i}$）を含めて\n",
    "\n",
    "$$y_{it}=\\alpha+\\beta x_{it}+\\gamma\\bar{x}_i+\\delta z_{i} +v_{it}\\qquad\\text{(式６)}$$\n",
    "\n",
    "をRE推定しても次の結果は成立する。\n",
    "\n",
    "$$\\hat{\\beta}_{FE}=\\hat{\\beta}_{CRE}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "この結果を利用して，以下では次の変数を加えて回帰式を設定する。\n",
    "* 時間に対して一定な変数\n",
    "    * `educ`, `black`，`hisp`\n",
    "* 一定ではないが，FEにいれると推定できなかった変数\n",
    "    * `exper`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "formula_cre2 = 'lwage ~ 1 + married + union + expersq \\\n",
    "                          + exper + educ + black + hisp \\\n",
    "                          + married_mean + union_mean + expersq_mean \\\n",
    "                          + d81 + d82 + d83 + d84 + d85 + d86 + d87'\n",
    "\n",
    "result_cre2 = RandomEffects.from_formula(formula_cre2, data=wagepan).fit()\n",
    "\n",
    "print(result_cre2.summary.tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "もう一度，FE対RE検定を行ってみよう。検定統計量（Statistics）は減少したが，以前として帰無仮説は高い優位水準で棄却できる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "result_cre2.wald_test(formula=restriction)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "## モデルの比較"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "パネル・データを扱う場合の通常のアプローチは，使える推定法を使いその結果を比べることから始める。以下では以下のモデルを比べる。\n",
    "* 通常のOLS\n",
    "* 固定効果モデル\n",
    "* ランダム効果モデル\n",
    "* 相関ランダムモデル"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true,
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "### OLS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`linearmodels`のモジュール`PooledOLS`では，観察単位や時間の区別なく全てのデータをプールしてOLS推定する。これは通常のOLSと同じ推定法と等しい。`PooledOLS`の関数`from_formula`を使い，以下のように推定する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "formula_pool = 'lwage ~ 1 + married + union + expersq \\\n",
    "                        + exper + educ + black + hisp \\\n",
    "                        +d81+d82+d83+d84+d85+d86+d87'\n",
    "\n",
    "result_pool = PooledOLS.from_formula(formula_pool, data=wagepan).fit()\n",
    "\n",
    "print(result_pool.summary.tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true,
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "### 比較表の作成"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "今までの推定結果を表にまとめるために，`linearmodels`の関数`compare`を使う。\n",
    "\n",
    "表作成の順番：\n",
    "1. `key`が表示したい推定方法の名前，`value`がその上で得た推定結果となる辞書を作る\n",
    "1. その辞書を`compare`の引数としてつかう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "res = {'Pooled OLS':result_pool,\n",
    "       'FE': result_fe,\n",
    "       'RE': result_re,\n",
    "       'CRE': result_cre2\n",
    "      }\n",
    "\n",
    "# compare(res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "このままでも良いか，この方法では推定結果を表示する順番を指定できない。例えば，`OLS`，`FE`，`RE`，`CRE`の順番で左から並べたいとしよう。その場合，`collections`パッケージにある関数`OrderedDict`をつかう。`{}`の中で並べた順番をそのまま維持してくれる関数である。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true,
    "tags": [
     "output_scroll"
    ]
   },
   "outputs": [],
   "source": [
    "res_ordered = OrderedDict(res)\n",
    "print(compare(res_ordered))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "* `married`\n",
    "    * OLSでの結婚プレミアムは10％以上あるが，FEでは半減している。これは観察単位の異質性$a_i$に生産性が含まれており，「生産性が高い人（高い$a_i$）は，結婚している可能性が高い」という考えと一貫性がある。即ち，$\\text{Cov}\\left(a_i,x_{it}\\right)>0$となり，これによりOLSでは以下のようなバイアスが発生すると解釈できる。\n",
    "        * 既婚 $\\Rightarrow$ 生産性（$a_i$）が高い $\\Rightarrow$ 賃金が上振れしやすい\n",
    "        * 未婚 $\\Rightarrow$ 生産性（$a_i$）が低い $\\Rightarrow$ 賃金が下振れしやすい\n",
    "    * FE推定法が示しているのは，そのようなバイアスを取り除いても，結婚プレミアムは存在する。考えられる理由は：\n",
    "        * 結婚は生産性を上昇させる\n",
    "        * 結婚は安定した生活を意味し，それに対して企業はより高い賃金を払う\n",
    "    * REは推定値はFEに比較的に近い。これは$\\hat{\\theta}>0.5$の値にも現れている。\n",
    "* `union`\n",
    "    * FEの値は労働組合の賃金に対する影響力を示している。\n",
    "    * OLSとFEの推定値を比較すると，約0.1減少している。OLSは労働組合の影響力を過大評価しており，観察単位の異質性が大きく働いていることがわかる。\n",
    "* `educ`，`black`，`hisp`\n",
    "    * OLSもREも推定値は似ている。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true,
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "### FD vs. FE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "1階差分モデルと固定効果モデルを比較する。\n",
    "* $T=2$の場合，FDとFEは同じ（if there is intercept in FE）\n",
    "* $T\\geq 3$の場合：\n",
    "    * GM仮定に対応する仮定の下ではFDもFEも不偏性・一致性を満たす。\n",
    "    * 誤差項に系列相関がない場合，FEの方が効率性が高い\n",
    "    * 誤差項の系列相関がある場合，FDの方が良い。\n",
    "        * 例えば，誤差項がランダム・ウォークの場合$\\Delta u_{it}$は系列相関はない。\n",
    "        * FD推定をして$\\Delta u_{it}$を検定する。\n",
    "    * $N$が少なく$T$が大きい場合（例：$N=20$と$T=30$），時系列の特性が強くなるので，FDの方が良い\n",
    "* 実証研究では，FDとFEの結果の両方を報告すること。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true,
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "### FE vs. RE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "固定効果モデルとランダム効果モデルを比較する。\n",
    "1. $a_i$はランダムか？\n",
    "    * 経済学の場合，説明変数は何らかの選択の結果の場合が多い。さらに，その選択が観察単位の特徴に依存する場合，$\\text{Cov}\\left(a_ix_{it}\\right)\\neq 0$となり，FEモデルの方が適切となる。\n",
    "    * 都道府県データのような場合，「大きな」母集団からランダムに抽出された説明変数とはならないかも知れない。むしろ，都道府県の切片がランダムではなく，単に異なると仮定する方が自然かも知れない。\n",
    "1. Hausman検定であれ他の検定であれ，間違う確率は存在する。以下では間違った場合どうなるかを考えた。\n",
    "    * $\\text{Cov}\\left(a_ix_{it}\\right)= 0$，即ち，REモデルが正しい場合：\n",
    "        * 誤差項の自己相関があるが，FE推定量は不偏性を満たす\n",
    "    * $\\text{Cov}\\left(a_ix_{it}\\right)\\neq 0$，即ち，FEモデルが正しい場合：\n",
    "        * GM仮定４が満たされないため，RE推定量は不偏性を満たさない\n",
    "\n",
    "Pros for RE\n",
    "1. FEの場合，時間に対して変化しない変数の係数を推定できない\n",
    "2. 回帰式が非線形の場合（例えば，Probit），FEでは対応できない。\n",
    "\n",
    "\n",
    "（結論）一般的にFEの方が適切な場合が多いのではないか。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "## シミュレーション"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "（目的）真のモデルでは観察単位の固定効果がある場合を考え，FE推定量とRE推定量を比較する。\n",
    "\n",
    "単回帰を想定する。以下を真の値として設定する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "b0 = 1  # 定数項\n",
    "b1 = 2  # スロープ係数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 推定値の計算"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "シミュレーションの基本パラメータ等の設定"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "N = 100  # 観察単位の数\n",
    "T = 5  # 時間数\n",
    "ai = np.linspace(0,10,N)  # 観察単位の異質性"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`for`ループによる`DataFrame`の作成。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "df_sim = pd.DataFrame()  # 空のDataFrame\n",
    "\n",
    "for (idx,a) in enumerate(ai):\n",
    "    \n",
    "    x = norm.rvs(a,1,size=T)  # T個のaが平均となるランダムな数\n",
    "    u = norm.rvs(0,1,size=T)  # 誤差項\n",
    "    y = b0 + b1*x + a + u  # 被説明変数\n",
    "    df_idx = pd.DataFrame({'id':[idx]*T,    # 観察個体のID\n",
    "                         'time':np.array(range(T))+2000,\n",
    "                         'Y':y,\n",
    "                         'X':x,\n",
    "                         'ai':[a]*T})\n",
    "    df_sim = pd.concat([df_sim,df_idx])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`DataFrame`の微調整。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# id と time の列を 整数型に変換（省いても問題ない）\n",
    "df_sim['id'] = df_sim['id'].astype(int)\n",
    "df_sim['time'] = df_sim['time'].astype(int)\n",
    "\n",
    "# MultiIndex化\n",
    "df_sim = df_sim.set_index(['id','time'])\n",
    "\n",
    "df_sim.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "固定効果モデルによる推定。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "form_sim_fe = 'Y ~ X + EntityEffects'\n",
    "\n",
    "sim_fe = PanelOLS.from_formula(form_sim_fe, data=df_sim).fit()\n",
    "\n",
    "print(sim_fe.summary.tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "ランダム効果モデルによる推定。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "form_sim_re = 'Y ~ 1 + X'\n",
    "\n",
    "sim_re = RandomEffects.from_formula(form_sim_re, data=df_sim).fit()\n",
    "\n",
    "print(sim_re.summary.tables[1])\n",
    "print('theta:', sim_re.theta.iloc[0,:].values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "$\\text{Cov}\\left(a_i,x_{it}\\right)>0$により上方バイアスが発生している。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "相関ランダム効果モデルによる推定。\n",
    "\n",
    "まず，観察単位の`X`の平均の列を追加する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "df_sim = add_col_mean(df_sim, 'X', 'X_mean')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "form_sim_cre = 'Y ~ 1 + X + X_mean'\n",
    "\n",
    "sim_cre = RandomEffects.from_formula(form_sim_cre, data=df_sim).fit()\n",
    "\n",
    "print(sim_cre.summary.tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "通常のOLS（Pooled OLS）による推定はバイアスが発生する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "form_sim_pool = 'Y ~ 1 + X'\n",
    "\n",
    "sim_pool = PooledOLS.from_formula(form_sim_pool, data=df_sim).fit()\n",
    "\n",
    "print(sim_pool.summary.tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 推定値の分布"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "基本的に上のシミュレーションのコードを応用する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "N = 100  # 観察単位数\n",
    "T = 5  # 年数\n",
    "R = 100  # シミュレーションの回数\n",
    "ai = np.linspace(0,10,N)  # 観察単位の異質性"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "（下のコードについて）\n",
    "* ランダム効果モデルの推定には`linearmodels`を使っている。固定効果モデルと通常のOLS推定にも`linearmodels`を使うと必要がない統計量も計算するため計算に時間が掛かる。少しでも計算時間を縮めるために「手計算」をする。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "bhat_fe_list = []  # FE推定値を入れるための空のリスト\n",
    "bhat_re_list = []  # RE推定値を入れるための空のリスト\n",
    "bhat_pool_list = []  # Pooled OLS推定値を入れるための空のリスト\n",
    "\n",
    "for _ in range(R):  # Rの値は下のコードで使わないので\"_\"に設定する\n",
    "    \n",
    "    df = pd.DataFrame()  # 空のDataFrame\n",
    "\n",
    "    # データの生成\n",
    "    for (idx,a) in enumerate(ai):\n",
    "        x = norm.rvs(a,1,size=T)  # T個のaが平均となるランダムな数\n",
    "        u = norm.rvs(0,1,size=T)  # T個の誤差項\n",
    "        y = b0 + b1*x + a + u  # T個の被説明変数\n",
    "        df_idx = pd.DataFrame({'id':[idx]*T,      # DataFrameへ表の追加\n",
    "                             'time':np.array(range(T))+2000,\n",
    "                             'Y':y,\n",
    "                             'X':x,\n",
    "                             'ai':[a]*T})\n",
    "        df = pd.concat([df,df_idx])\n",
    "        \n",
    "    # RE推定\n",
    "    df_re = df.set_index(['id','time'])  # MultiIndex化\n",
    "    form_sim_re = 'Y ~ 1 + X'\n",
    "    sim_re = RandomEffects.from_formula(form_sim_re, data=df_re).fit()\n",
    "    bhat_re_list.append(sim_re.params.iloc[1])\n",
    "\n",
    "    # FE推定\n",
    "    df_fe = df.loc[:,['Y','X']] - df.groupby('id')[['Y','X']].transform('mean')\n",
    "    Yfe = df_fe.loc[:,'Y'].values\n",
    "    Xfe = df_fe.loc[:,'X'].values[:,None]  # [:,None]は(N*T,1)の行列に変換\n",
    "    bhat_fe = (np.linalg.inv(Xfe.T@Xfe)@Xfe.T@Yfe)[0]\n",
    "    bhat_fe_list.append(bhat_fe)\n",
    "    \n",
    "    # Pooled OLS推定\n",
    "    c = np.ones(N*T)\n",
    "    Xpool = np.stack([c, df.loc[:,'X'].values], axis=1)\n",
    "    Ypool = df.loc[:,'Y'].values\n",
    "    bhat_pool = (np.linalg.inv(Xpool.T@Xpool)@Xpool.T@Ypool)[1]\n",
    "    bhat_pool_list.append(bhat_pool)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "分布の図示"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "xx=np.linspace(1.6,3.1,num=100)  # 図を作成するために横軸の値を設定\n",
    "\n",
    "kde_model_fe=gaussian_kde(bhat_fe_list)  # FE推定量のカーネル密度関数を計算\n",
    "\n",
    "kde_model_re=gaussian_kde(bhat_re_list)  # RE推定量のカーネル密度関数を計算\n",
    "\n",
    "kde_model_pool=gaussian_kde(bhat_pool_list)  # Pooled OLS推定量のカーネル密度関数を計算\n",
    "\n",
    "plt.plot(xx, kde_model_fe(xx), 'g-', label='FE')  # FE推定量の分布プロット\n",
    "plt.plot(xx, kde_model_re(xx),'r-', label='RE')  # RE推定量の分布プロット\n",
    "plt.plot(xx, kde_model_pool(xx),'k-', label='Pooled OLS')  # RE推定量の分布プロット\n",
    "plt.axvline(x=b1,linestyle='dashed')  # 真の値での垂直線\n",
    "plt.ylabel('Kernel Density')  # 縦軸のラベル\n",
    "plt.legend()  # 凡例\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "## 標準誤差の問題"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "(式３)の固定効果モデルを考えよう。パネル・データの場合，次の2つの問題が発生する場合がある。\n",
    "* 残差の不均一分散\n",
    "    * 固定効果モデルは，変数の変換後に通常のOLS推定をおこなうため，残差に不均一分散が存在するかも知れない。（[不均一分散](chap:14-hetero)を参照。）\n",
    "* 残差の相関\n",
    "    * 直観的に説明するために次の例を考えよう。賃金のパネル・データを作成する際に，まず複数の会社をランダムに抽出し，それらの会社から複数の労働者をランダム抽出することにより，各労働者の時系列データを収集するとしよう。パネル・データ分析をすると次の残差の相関が発生する可能性がある。\n",
    "    * **＜個体内での残差の相関＞** 労働者Aさんを考えよう。理想的な仮定の下では，Aさんの残差は時系列的に独立となる。しかし，その仮定が満たされない場合，Aさんの残差には自己相関が存在することになる。Bさんも，Cさんも，，，Zさんにも同様のことが言える。このタイプの相関は，**個体内**での期間をまたがる残差の相関である。\n",
    "    * **＜個体間での残差の相関＞** ある期間を考えよう。会社Firm0で働くA0さん，B0さん，C0さん，，，Z0さんは，同じ会社で働いているため共通の特徴があり，残差は相関するかも知れない。別の会社Firm1で働くA1さん，B1さん，C1さん，，，Z1さんも，同様に残差は相関する可能性がある。このタイプの相関は，期間内での**個体間**の残差の相関となる。\n",
    "\n",
    "これらの問題が発生しても，基本的な仮定の下（特に，[GM仮定４が満たされない場合](chap:18-GM4)で詳しく検討するガウス・マルコフ定理４が満たされる場合）でFE推定量は不偏性と一致性を満たす。しかし，係数の標準誤差は有効ではなくなり検定が無効となる。その対処方法として，`linearmodels`には主に次の２つの方法が用意されている。\n",
    "\n",
    "1. 不均一分散頑健推定\n",
    "    * 詳しくは，[不均一分散](chap:14-hetero)を参照。\n",
    "1. クラスター頑健推定\n",
    "    * クラスターとは何らかの特徴を共有する観測値がグループ化された集まりである。クラスター頑健推定では，（１）残差の不均一分散と（２）格クラスター内での残差の相関の両方を考慮した頑健推定となる。（２）には上で説明した次のパターンが挙げられる。\n",
    "        * 個体**内**での残差の自己相関であり，時系列的な相関に該当する。\n",
    "        * 個体**間**での残差の相関であり，クロス・セクションの相関に該当する。\n",
    "    \n",
    "以下では，この２つの頑健推定の使い方を説明する。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true,
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "### 分散の確認"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`linearmodels`には残差の均一分散を調べるBreusch-Pagan検定やWhite検定をおこなうメソッドは用意されていない。ここでは`statsmodels`を使いWooldridge (2016,p.253)で説明されているWhite検定に基づく検定をおこなう。次に必要な変数を作成する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# 残差\n",
    "u_hat = result_fe.resids.values.flatten()\n",
    "\n",
    "# 被説明変数の予測値\n",
    "y_hat = result_fe.fitted_values.values.flatten()\n",
    "\n",
    "# DataFrameの作成\n",
    "df_white = pd.DataFrame({'u_hat':u_hat,'y_hat':y_hat})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "（上のコードの説明）\n",
    "* `resids`は残差を`DataFrame`として取得する属性\n",
    "* `fitted_values`は予測値を`DataFrame`として取得する属性\n",
    "* `values`は`DataFrame`を`array`として返す属性\n",
    "* `flatten()`は`array`が2次元になっているのを1次元に変換するメソッド\n",
    "    * `array([[..],[..],...[...]])`を`array([....])`に変換する。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "検定に使う式\n",
    "\n",
    "$$\\hat{u}^2=\\beta_0+\\beta_1\\hat{y}+\\beta_2\\hat{y}^2+e$$\n",
    "\n",
    "* $\\text{H}_0：\\beta_1=\\beta_2=0$（均一分散）\n",
    "* $\\text{H}_A：$帰無仮説は成立しない"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "form_white = 'I(u_hat**2) ~ y_hat + I(y_hat**2)'\n",
    "\n",
    "res_white = smf.ols(form_white, data=df_white).fit()\n",
    "\n",
    "print(res_white.summary().tables[1])\n",
    "print('F検定のp値：',res_white.f_pvalue)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "帰無仮説は棄却できない。\n",
    "\n",
    "次に図を使い確認する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "b0 = res_white.params.iloc[0] # beta0\n",
    "b1 = res_white.params.iloc[1] # beta1\n",
    "b2 = res_white.params.iloc[2] # beta2\n",
    "\n",
    "xx = np.linspace(min(y_hat), max(y_hat), 100) # x軸の値\n",
    "z = b0 + b1*xx + b1*xx**2  # 検定に使った式\n",
    "\n",
    "plt.scatter(y_hat, u_hat**2)  # u_hat, y_hatの散布図\n",
    "plt.plot(xx, z, 'red', linewidth=3)  # 検定の式の曲線\n",
    "plt.xlabel('y_hat')\n",
    "plt.ylabel('u_hat^2')\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "（解説）\n",
    "* 上の検定で推定した式は赤い線である。殆ど平行になっているため帰無仮説を棄却できなかった。\n",
    "* 図の中で`y_hat`の値が`-0.2`を境に`u_hat`の変動はより大きくなっており，不均一分散の疑いが高い。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "```{note}\n",
    "残差の自己相関の検定については後日追加する予定。\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true,
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "### 対処方法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true,
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "#### 説明と計算"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "＜不均一分散頑健的推定の使い方＞\n",
    "* メソッド`fit()`に以下の引数を指定する。\n",
    "    * クラスター頑健推定を指定：`cov_type='robust'`\n",
    "\n",
    "＜クラスター頑健的推定の使い方＞\n",
    "* メソッド`fit()`に以下の引数を指定する。\n",
    "    * クラスター頑健推定を指定：`cov_type='clustered'`\n",
    "    * グループ化するクラスターを指定：\n",
    "        * 個体のクラスターを指定する場合：`cluster_entity=True`\n",
    "        * 時間のクラスターを指定する場合：`cluster_time=True`\n",
    "        * 個体と時間の２次元のクラスターを指定する場合：`cluster_entity=True`及び`cluster_time=True`\n",
    "    * クラスター数が小さい場合の影響を調整する：`group_debias=True`（デフォルトは`False`だが，`True`を推奨する）\n",
    "\n",
    "（注意）\n",
    "* 係数の推定値は変わらない。\n",
    "* 係数の標準誤差だけが修正され，有効な検定ができるようになる。\n",
    "\n",
    "\n",
    "`wagepan`を使って上で推定した式にクラスター頑健的推定を使う。\n",
    "\n",
    "＜不均一分散頑健的推定の場合＞"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "mod_fe_clus1 = PanelOLS.from_formula(formula_fe, data=wagepan)\n",
    "\n",
    "res_fe_clus1 = mod_fe_clus1.fit(cov_type='robust')\n",
    "\n",
    "print(res_fe_clus1.summary.tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "＜クラスター頑健推定の場合＞"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "mod_fe_clus1 = PanelOLS.from_formula(formula_fe, data=wagepan)\n",
    "\n",
    "res_fe_clus1 = mod_fe_clus1.fit(cov_type='clustered',\n",
    "                                cluster_entity=True,\n",
    "                                group_debias=True)\n",
    "\n",
    "print(res_fe_clus1.summary.tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "不均一分散頑健推定とクラスター頑健推定の係数の標準誤差（従って，`t`値）は異なる。理由は，残差にある個体ごとの自己相関の存在である。実際，残差の自己相関AR(1)を検定すると，その存在を棄却できない。従って，不均一分散頑健推定は自己相関を取りこぼすことになり，標準誤差が小さくなっている。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "````{note}\n",
    "固定効果モデルを使いクラスター頑健推定について説明したが，ランダム効果モデルや通常の回帰分析でも使うことができる。例えば，通常の回帰分析の場合は次のコードとなる。\n",
    "```\n",
    "formula = 'lwage ~ 1 + married + union + expersq \\\n",
    "                   +d81+d82+d83+d84+d85+d86+d87'\n",
    "mod = PooledlOLS.from_formula(formula, data=wagepan)\n",
    "res = mod.fit(cov_type='clustered',\n",
    "              cluster_entity=True,\n",
    "              group_debias=True)\n",
    "```\n",
    "クラスター頑健推定の妥当性は，データの収集プロセスなどを考慮して検討する必要があるだろう。\n",
    "````"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "```{warning}\n",
    "クラスター頑健推定について次の点に注意しよう。\n",
    "* 大標本（$N\\rightarrow\\infty$）の下で有効となる。\n",
    "* 個体数が十分に大きい必要がある。個体数が小さい場合に使うと，問題を悪化させる可能性がある。\n",
    "* 個体数よりも時系列のデータが長いパネル・データの場合は意味がない推定となってしまう。\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "#### `p`値について"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "一方で，`linearmodels`でクラスター頑健推定を使う場合，`p`値の計算については以下で説明する点について覚えておくと良いだろう。次の変数を定義しよう。\n",
    "\n",
    "* $n$：各個体の観測値\n",
    "* $G$： 個体数（クラスターの数，Groupの**G**）\n",
    "* $N=nG$: 標本の大きさ\n",
    "* $k$: 定数項以外の説明変数の数\n",
    "\n",
    "`linearmodels`は，クラスター頑健推定の`p`値を計算する際，`t`分布に自由度`N-G-k`を使っている。変数`married`を使って，数値で確認してみる。`N-G-k`は推定結果のオブジェクト`res_fe_clus1`の属性`.df_resid`かる抽出することができる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# 自動計算のp値\n",
    "pval_auto_calc = res_fe_clus1.pvalues['married']\n",
    "\n",
    "# 「手計算」のp値\n",
    "dof = res_fe_clus1.df_resid  # 自由度\n",
    "pval_manual_calc = 2*(\n",
    "    1-t.cdf(res_fe_clus1.tstats['married'], dof) )\n",
    "\n",
    "print('marreidのp値\\n'+'-'*27)\n",
    "print(f'自動計算：　　{pval_auto_calc}')\n",
    "print(f'「手動計算」：{pval_manual_calc}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "同じ値だということが分かる。\n",
    "\n",
    "一方，Stata（`xtreg`コマンド）では，この場合の自由度は`G-1`を使っている。この値を使って`p`値を計算し直して比較してみよう。推定結果のオブジェクト`res_fe_clus1`の属性`entity_info`は，個体に関する情報を`Series`として返す。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "res_fe_clus1.entity_info"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "データセットはbalanced dataset（$n_1=n_2=\\cdots=n_G$）であり，それぞれの個体には`8`の観測値があることが分かる。また`total`は個体の総数を表しているので，次の変数を作成しよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "G = int( res_fe_clus1.entity_info['total'] )\n",
    "\n",
    "# Stataで使う自由度\n",
    "dof_stata = G-1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "この値を使い`p`値を計算してみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "2 * ( 1-t.cdf(res_fe_clus1.tstats['married'], dof_stata) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "小数点第3位までは`linearmodels`と同じ値だが，それ以降はStataの結果と異なる。Stataの結果とは，若干異なる事を覚えておくと，後になって戸惑うことを避けることができるだろう。\n",
    "\n",
    "では，どちらを使えば良いのだろうか。Stataが使う自由度$G-1$の方が一般的であるため，その慣例に従えば良いだろう。一方で，自由度が`30`以上の`t`分布の場合，標準正規標分布との誤差は非常に小さな値になる。両方の自由度を表示してみる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "print(f'linearmodelsの自由度：{dof}')\n",
    "print(f'Stataの自由度：{dof_stata}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "２つのケースとも，自由度は大きいため，上の計算が示すように，`p`値の差は大きくない。実際に，標準正規分布を使った場合の`p`値を計算してみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "2 * ( 1-norm.cdf(res_fe_clus1.tstats['married']) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "近い値になっていることが分かる。従って，手元のデータが$N\\rightarrow\\infty$と近似でき，且つ，個体数$G$が十分に大きいと判断される場合は，次のように使い分けるのも一案だろう。\n",
    "* 暫定的な推定には，`linearmodels`のデフォルトの`p`値を使い，最終盤の推定式を模索する。\n",
    "* 卒論などの提出版には，自由度$G-1$とする`p`値を「手計算」する。\n",
    "\n",
    "いずれ`linearmodels`も変更するかもしれないが，現状では，この方法が最も現実的だろう。"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.9"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}