19_IV2SLS.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 操作変数法と２段階OLS"
   ]
  },
  {
   "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 wooldridge\n",
    "\n",
    "from linearmodels.iv import IV2SLS\n",
    "from scipy.stats import chi2, gaussian_kde, multivariate_normal, norm, uniform\n",
    "\n",
    "# 警告メッセージを非表示\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**＜仮定４が満たされない場合＞**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## 説明"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 基本的な考え方"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "仮定４：Zero Conditional Mean $\\text{E}(u|x)=0$\n",
    "\n",
    "仮定４a：$\\text{Cov}(x,u)=0$\n",
    "\n",
    "この仮定が満たされない理由に次の３つがあり（他にもある），その対処法となる推定方法について説明する。\n",
    "* 欠落変数\n",
    "* 同時方程式\n",
    "* 測定誤差\n",
    "\n",
    "この仮定４もしくは４aが満たされている場合，$x$は**外生変数**という。一方，この仮定が満たされない場合，OLS推定量は不偏性も一致性も満たさない。その場合の説明変数は**内生変数**とよばれる。説明変数の内生性は経済問題に多く存在すると考えられる。例えば，既出の賃金関数を考えよう。賃金は教育や経験などに依存しているが，労働者の内在的な能力にも依存していると考えられる。能力を捉える変数が回帰式にない場合（欠落変数），その効果は誤差項$u$に入ることになる。もちろん，能力が労働者の中でランダムに存在し，仮定４もしくは４aを満たしているのであれば問題がない。しかし説明変数である教育が誤差項に含まれる能力と何らかの関係がある場合（例えば，能力の高い人がより高い教育水準を選択する），仮定４もしくは４aは満たされないことになり，OLS推定量は不偏性を満たさない（**欠落変数バイアス**がある）。\n",
    "\n",
    "```{note}\n",
    "内生変数とは誤差項と相関する変数を意味する。従って，被説明変数も内生変数である。しかし，以下では「内生変数」という用語を「内生的説明変数」の省略形として使うことにする。\n",
    "```\n",
    "\n",
    "このような場合に役に立つのが**操作変数法**（`I`nstrumental `V`ariable Estimation）と呼ばれる手法である。この推定法では，ある条件を満たす内生的説明変数の代わりになる**操作変数**（外生変数）を使うことにより，一致性を満たす推定量を得ることが可能となる。\n",
    "\n",
    "（注意点）\n",
    "* IV推定量は**一致性**を満たすが，この特性を活かすためには十分に大きな標本が必要である。\n",
    "* 標本の大きさが小さい場合，IV推定量は不偏性を失う。\n",
    "* OLS推定量と比べてIV推定量の標準誤差は大きくなる（効率性が低い）。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "基本的なアイデアを整理するために次の単回帰式を考えよう。\n",
    "\n",
    "$$\n",
    "y=\\beta_0+\\beta_1w+u\n",
    "$$ (eq:19-1)\n",
    "\n",
    "* $w$は説明変数\n",
    "* $\\text{Cov}(w,u)\\neq 0$（即ち，$w$は内生変数）\n",
    "\n",
    "**＜操作変数の３つの条件＞**\n",
    "\n",
    "更に，$w$に以下の条件を満たす**操作変数**（instruments）$z$があるとしよう。\n",
    "1. $z$は回帰式に含まれない（除外条件）\n",
    "1. $\\text{Cov}(z,w)\\neq 0$（妥当性; $w$と高い相関関係がある）\n",
    "1. $\\text{Cov}(z,u)=0$（外生性; 回帰式の誤差項と無相関）\n",
    "\n",
    "この場合，操作変数法を用いて（大標本のもとで）一致性を満たす$\\hat{\\beta}_1$を得ることが可能となる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "＜推定方法の考え方＞\n",
    "\n",
    "操作変数（IV）推定方は**２段階最小２乗法**の特殊なケースとして理解できる。\n",
    "* 操作変数の数と内生変数の数が等しい場合の推定法を操作変数（IV）推定法\n",
    "* 操作変数の数が内生変数の数を上回る場合の推定法を２段階最小二乗法（2 Stage Least Squares; 2SLS）\n",
    "\n",
    "**第１段階OLS**\n",
    "* 次式をOLS推定する。\n",
    "\n",
    "$$w=\\pi_0+\\pi_1z+v$$\n",
    "\n",
    "* これにより$w$を予測値$\\hat{w}$と残差$\\hat{v}$に分解する。\n",
    "\n",
    "    $$\n",
    "    w=\\hat{w}+\\hat{v},\\qquad\\hat{w}=\\hat{\\pi}_0+\\hat{\\pi}_1z\n",
    "    $$ (eq:19-2)\n",
    "\n",
    "    1. $\\hat{w}$（予測値）：外生変数で説明できる$w$の部分（$z$と相関する部分）\n",
    "    1. $\\hat{v}$（残差）：残り全て（$u$と相関する$w$の部分は$\\hat{v}$に吸収される）\n",
    "    \n",
    "* $\\hat{\\pi}_1$の統計的優位性の確認\n",
    "    * 一般的に有効な操作変数は以下を満たす\n",
    "        * 操作変数の$t$値の絶対値 $>3.2$\n",
    "        * この目安を満たさなければ弱操作変数の可能性があり，その場合，推定量は不偏性・一致性を満たさない\n",
    "    * 帰無仮説$\\text{H}_0:\\;\\hat{\\pi}_1=0$が棄却され，上の基準をクリアすれば次のステップへ\n",
    "\n",
    "**第２段階OLS**\n",
    "* 予測値$\\hat{w}$を使い次式をOLS推定する。\n",
    "\n",
    "    $$\n",
    "    y=\\gamma_0+\\gamma_1\\hat{w}+e\n",
    "    $$ (eq:19-3)\n",
    "    \n",
    "* $\\hat{\\gamma}_1$：IV（instrumental variable）推定量"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "（注意点）\n",
    "\n",
    "* 「手計算」で第１・２段階を別々にOLS推定すると，$\\hat{\\gamma}_j$を得ることができるが，推定量の標準誤差，検定統計量，決定係数$R^2$は有効ではない。\n",
    "* `Python`の専用パッケージを使と，\n",
    "    * 大標本のもとで推定量と$t$値は有効となる。\n",
    "    * $R^2$は特に有用な情報を提供しない（**マイナス**になり得る）。\n",
    "    * パラメータ制約を検定する$F$検定をする場合は「手計算」ではなくパッケージで提供されたコマンドを使うこと。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### `linearmodels`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "**＜`linearmodels`の使い方＞**\n",
    "\n",
    "回帰式を文字列で表し操作変数法を使うためにはサブパッケージ`IV2SLS`の`from_formula`をインポートする。回帰式の文字列に関して上述した定数項についての違い以外は`statsmodels`と同じである。ただし，下の回帰式の一般型に沿って内生変数と操作変数を`~`で挟んで`[ ]`の中に入れる。\n",
    "\n",
    "**＜`IV2SLS`回帰式の一般形＞**\n",
    "```\n",
    "被説明変数 ~ 定数項 + 外生変数 + [内生変数 ~ 操作変数]\n",
    "```\n",
    "* 定数項，外生変数がない場合は省いても良い。\n",
    "* `[ ]`がない場合は通常のOLSとして計算される。その際，以前説明したメソッド`fit()`のオプションに要注意。\n",
    "* 外生変数，操作変数は複数でも可"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ケース１：単純な`IV`推定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(sec:19-case1)=\n",
    "## ケース１：単純な`IV`推定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### データ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "既婚女性の教育の収益率に関するデータ`mroz`を利用して使い方を説明する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true,
    "tags": [
     "output_scroll"
    ]
   },
   "outputs": [],
   "source": [
    "mroz = wooldridge.data('mroz').dropna(subset=['lwage'])  # 列'lwage'にNaNがある行は削除する\n",
    "\n",
    "wooldridge.data('mroz',description=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "ケース１では以下の場合を考える。\n",
    "* 被説明変数：`lwage`（既婚女性の賃金; 対数）\n",
    "* 内生的説明変数：`educ`（既婚女性の教育年数）\n",
    "* 操作変数：`fatheduc`（父親の教育年数）\n",
    "* 外生的説明変数：なし\n",
    "\n",
    "（考え方）\n",
    "\n",
    "誤差項に既婚女性の能力が含まれている可能性があるため`educ`は内生変数の疑いがある。父親の教育年数`fatheduc`は既婚女性の教育年数`educ`と正の相関性があると思われる一方，能力自体とは無相関と仮定。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### 自動計算"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "IV推定法は2SLSの特殊なケースとして２ステップで推定することを説明したが，ここでは自動的に２ステップを計算する場合を紹介する。\n",
    "\n",
    "まず回帰式を決める。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "form_1 = 'lwage ~ 1 + [educ ~ fatheduc]'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`IV2SLS`モジュールの`from_formula`を使うことにより，`statsmodels`の`ols`のように回帰式を文字列で指定できる。次式では推定するモデルを設定する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "mod_1 = IV2SLS.from_formula(form_1, data=mroz)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`statsmodels`の`ols`のようにメソッド`.fit()`を使い推定する。（以前説明したオプションについての説明を参照）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "res_1 = mod_1.fit(cov_type='unadjusted')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`res_1`の属性`.summary`，さらに`summary`の属性`tables`を使ってパラメータの部分だけを表示する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "print(res_1.summary.tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`educ`の推定値は，上で説明したように`fatheduc`を操作変数として計算した結果である。$p$値は`0.1`よりも低いため`10%`有意水準では係数`0`の帰無仮説を棄却できるが，`5%`水準では棄却できない。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### 「手計算」"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "第１・２段階のOLS推定を別々に試みる。\n",
    "\n",
    "---\n",
    "**第１段階のOLS**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "stage_1 = 'educ ~ 1 + fatheduc'  # 回帰式\n",
    "\n",
    "res_stage_1 =IV2SLS.from_formula(stage_1, data=mroz).fit(cov_type='unadjusted')  # OLS推定\n",
    "\n",
    "mroz['educ_fit'] = res_stage_1.fitted_values  # educの予測値を取得"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "上の３行目のでは`res_stage_1`の属性`.fitted_values`を使い予測値を取得している。`statsmodels`を使いOLS推定した際に使った`.fittedvalues`と異なるメソッド名になっていることに注意しよう。\n",
    "\n",
    "---\n",
    "**第２段階のOLS**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "stage_2 = 'lwage ~ 1 + educ_fit'\n",
    "\n",
    "res_stage_2 =IV2SLS.from_formula(stage_2, data=mroz).fit(cov_type='unadjusted')  # OLS推定\n",
    "\n",
    "print(res_stage_2.summary.tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "「自動計算」の場合と比べると，`Parameter`は同じことが確認できる。しかし`Std. Err.`は異なり，それに基づく他の推定値も異なることに注意。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### $OLS$推定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "確認のために，操作変数法を使わずに直接OLS推定をおこなうとどうなるかを確認しよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "form_ols = 'lwage ~ 1 + educ'\n",
    "\n",
    "res_ols =IV2SLS.from_formula(form_ols, data=mroz).fit()  # OLS推定\n",
    "\n",
    "print(res_ols.summary.tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "操作変数法の推定量と大きく異なることが分かる。操作変数法と比べて`educ`の$OLS$推定量は約２倍になり，既婚女性の教育の収益率を過大評価している。`educ`と誤差項に相関性がると推測できる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## ケース２：より複雑な`IV`推定（2SLS）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 説明"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "ケース１を次のように拡張する。\n",
    "* 複数の外生変数の導入（`exper`，`expersq`）\n",
    "* 複数の操作変数の導入（`fatheduc`，`motheduc`）\n",
    "\n",
    "このように複雑化しても基本的な考え方は同じである。\n",
    "\n",
    "$$\n",
    "y=\\beta_0+\\beta_1x_1+\\beta_2x_2+\\beta_3w+u\n",
    "$$ (eq:19-4)\n",
    "    \n",
    "* $\\text{Cov}(w,u)\\neq 0$（$w$は内生変数）\n",
    "* $\\text{Cov}(x_k,u)=0,\\;k=1,2$（$x_k$は外生変数）\n",
    "\n",
    "$w$の操作変数$z_j,\\;j=1,2$は次の条件を満たす必要がある。\n",
    "* $\\text{Cov}(z_j,w)\\neq 0$\n",
    "* $\\text{Cov}(z_j,u)=0$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "＜2SLSの考え方＞\n",
    "\n",
    "**第１段階OLS**\n",
    "* 次式をOLS推定する。\n",
    "\n",
    "    $$\n",
    "    w=\\pi_0+\\pi_1z_1+\\pi_2z_2+\\pi_3x_1+\\pi_4x_2+v\n",
    "    $$ (eq:19-5)\n",
    "    \n",
    "    * $x_1$と$x_2$が入る直感的な理由\n",
    "        * 下で説明するように，式[](eq:19-4)の`w`の代わりに式[](eq:19-4)を使って計算する`w`の予測値を代わりに使うことになる。その際，式[](eq:19-4)に$x_1$と$x_2$が入っていないと，`w`に元々あった$x_1$と$x_2$との関係性を除外することになる。式[](eq:19-4)に$x_1$と$x_2$を入れるのは$x_1$と$x_2$との相関を考慮した`w`の予測値にするためである。\n",
    "* $w$を２つの要素（$\\hat{w}$と$\\hat{v}$）に分解\n",
    "\n",
    "    $$\n",
    "    w=\\hat{w}+\\hat{v},\n",
    "    \\qquad\n",
    "    \\hat{w}=\\hat{\\pi}_0+\\hat{\\pi}_1z_1+\\hat{\\pi}_2z_2+\\hat{\\pi}_3x_1+\\hat{\\pi}_4x_2\n",
    "    $$ (eq:19-6)\n",
    "    \n",
    "    1. $\\hat{w}$（予測値）：外生変数だけで説明される$w$の部分（$x_1$，$x_2$，$z_1$，$z_2$と相関する部分）\n",
    "    1. $v$（残差）：残り全て（$u$と相関する$w$の部分はこれに吸収される）\n",
    "    \n",
    "* $\\hat{\\pi}_1$と$\\hat{\\pi}_2$の優位性の確認\n",
    "    * 操作変数の$F$値 $>10$\n",
    "    * この目安を満たさなければ弱操作変数の可能性\n",
    "        * 推定量は不偏性・一致性を満たさない\n",
    "    * $\\text{H}_0:\\;\\hat{\\pi}_1=\\hat{\\pi}_2=0$が棄却され，上の基準をクリアすれば次のステップへ\n",
    "\n",
    "**第２段階OLS**\n",
    "* 予測値$\\hat{w}$を使い次式をOLS推定する。\n",
    "\n",
    "    $$\n",
    "    y=\\gamma_0+\\gamma_1x_1+\\gamma_2x_2+\\gamma_3\\hat{w}+e\n",
    "    $$ (eq:19-7)\n",
    "    \n",
    "* IV（instrumental variable）推定量：$\\hat{\\gamma}_3$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### `educ`と操作変数の相関性チェック"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "内生的説明変数と操作変数のOLS推定を使い，相関性の検定をおこなう。\n",
    "\n",
    "上述のとおり，一般的に有効な操作変数は以下を満たす。\n",
    "* 操作変数が１つの場合\n",
    "    * 操作変数の$t$値の絶対値 $>3.2$\n",
    "* 複数の操作変数の場合\n",
    "    * 操作変数の$F$値 $>10$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "form_check = 'educ ~ 1 + fatheduc + motheduc'\n",
    "\n",
    "res_check =IV2SLS.from_formula(form_check, data=mroz).fit(cov_type='unadjusted')\n",
    "\n",
    "print(res_check.summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "* `fatheduc`と`motheduc`のそれぞれの係数\n",
    "    * `T-stat`（$t$値）$>3.2$\n",
    "    * `P-value`（$p$値）は約$0$\n",
    "* ２つの操作変数の係数が同時に$0$という帰無仮説の検定\n",
    "    * `F-statistic`（$F$値）$>10$\n",
    "    * `P-value(F-stat)`（$t$値）も約$0$\n",
    "\n",
    "従って，`educ`と操作変数の相関性は高い。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### IV推定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "上述した`IV2SLS`回帰式の一般形に基づいて回帰式を設定する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "form_2 = 'lwage ~ 1 + exper + expersq +[educ ~ fatheduc + motheduc]'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "操作変数法を使い推定"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "res_2 =IV2SLS.from_formula(form_2, data=mroz).fit(cov_type='unadjusted')\n",
    "\n",
    "print(res_2.summary.tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "IVが１つのケースと比べて`educ`の係数自体は大きく変わってはいないが，`5%`有意水準でも係数`0`の帰無仮説を棄却できるようになっている。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## ケース３：複数の内生変数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 説明"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "複数の内生変数がある場合を考える。まず重要なのは，**各内生変数に１つ以上の操作変数が存在**していることが前提となる。\n",
    "\n",
    "以下では次の仮定の下でシミュレーションを通して推定方法を説明する。\n",
    "* 母集団回帰式\n",
    "\n",
    "    $$\n",
    "    y=\\beta_0+\\beta_1x_1+\\beta_2x_2+\\beta_3x_3+\\gamma_1w_1+\\gamma_2w_2+e\n",
    "    $$ (eq:19-case3)\n",
    "    \n",
    "    * 誤差項$e$は全ての説明変数と無相関\n",
    "    \n",
    "* $w_1$は$x_1$と，$w_2$は$x_2$とそれぞれ相関する。\n",
    "    * $\\text{Cov}(w_1,x_1)=0.5$\n",
    "    * $\\text{Cov}(w_2,x_2)=0.5$\n",
    "    * 説明変数の他の共分散は全て`0`とする。\n",
    "    \n",
    "* $x_1$と$x_2$は欠落変数として次の回帰式を推定する。\n",
    "\n",
    "    $$\n",
    "    y_i=\\beta_0+\\beta_3x_{i3}+\\gamma_1w_{i1}+\\gamma_2w_{i2}+u_i\n",
    "    $$ (eq:19-case3a)\n",
    "    \n",
    "    * $u_i=e+\\beta_1x_{i1}+\\beta_2x_{i2}$\n",
    "    * $x_{i3}$は外生変数\n",
    "    * $w_{i1}$は内生変数，即ち，$\\text{Cov}(w_{i1},u_i)=\\text{Cov}(w_{i1},x_{i1})\\neq0$\n",
    "    * $w_{i2}$は内生変数，即ち，$\\text{Cov}(w_{i2},u_i)=\\text{Cov}(w_{i2},x_{i2})\\neq0$\n",
    "\n",
    "* $w_{ij}$の操作変数は$z_{ij}$，$j=1,2$とする。\n",
    "    * $\\text{Cov}(z_{ij},w_{ij})=0.8\\neq 0$（妥当性）\n",
    "    * $\\text{Cov}(z_{ij},u_i)=0$（外生性）\n",
    "\n",
    "* ここでガウス・マルコフ定理の仮定４が満たされない理由は，欠落変数$x_1$と$x_2$の存在である。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "内生変数が複数あっても推定の基本的な考え方は変わらない。違いは，第１段階において内生変数の数だけ（この例では２回）OLS推定をする必要がある点のみである。\n",
    "\n",
    "**第１段階OLS**\n",
    "* 次の２式をOLS推定する。\n",
    "\n",
    "    $$\n",
    "    \\begin{align*}\n",
    "    w_1&=a_{0}+a_{1}z_1+a_{2}z_2+a_{3}x_3+v_1\\\\\n",
    "    w_2&=b_{0}+b_{1}z_1+b_{2}z_2+b_{3}x_3+v_2\n",
    "    \\end{align*}\n",
    "    $$\n",
    "\n",
    "* それぞれを予測値$\\hat{w}_j$と残差$\\hat{v}_j$に分解する。\n",
    "\n",
    "    $$\n",
    "    \\begin{align*}\n",
    "    w_j&=\\hat{w}_j+\\hat{v}_j,\\quad j=1,2\\\\\n",
    "    \\hat{w}_1&=\\hat{a}_{0}+\\hat{a}_{1}z_1+\\hat{a}_{2}z_2+\\hat{a}_{3}x_3\\\\\n",
    "    \\hat{w}_2&=\\hat{b}_{0}+\\hat{b}_{1}z_1+\\hat{b}_{2}z_2+\\hat{b}_{3}x_3\n",
    "    \\end{align*}\n",
    "    $$\n",
    "\n",
    "    1. $\\hat{w}_j$（予測値）：外生変数で説明できる$w_j$の部分\n",
    "    1. $\\hat{v}_j$（残差）：残り全て（$u$と相関する$w$の部分は$\\hat{v}_j$に吸収される）\n",
    "    \n",
    "* $\\hat{a}_{1}$と$\\hat{b}_{2}$の統計的優位性の確認\n",
    "    * 一般的に有効な操作変数は以下を満たす\n",
    "        * 操作変数の$t$値の絶対値 $>3.2$\n",
    "        * この目安を満たさなければ弱操作変数の可能性があり，その場合，推定量は不偏性・一致性を満たさない\n",
    "    * 帰無仮説 $\\text{H}_0:\\;\\hat{a}_{1}=\\hat{b}_{2}=0$が棄却され，上の基準をクリアすれば次のステップへ\n",
    "\n",
    "**第２段階OLS**\n",
    "* 予測値$\\hat{w}_j$，$j=1,2$，を使い次式をOLS推定する。\n",
    "\n",
    "    $$\n",
    "    y=\\mu_0+\\mu_1\\hat{w}_1+\\mu_2\\hat{w}_2+\\mu_3x_3+\\varepsilon\n",
    "    $$\n",
    "    \n",
    "* $\\hat{\\mu}_1$と$\\hat{\\mu}_２$：IV（instrumental variable）推定量"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### データの作成"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "まずデータを作成するが，真のパラメターの値を次のように仮定する。\n",
    "\n",
    "$$\n",
    "\\beta_0=\\beta_1=\\beta_2=\\beta_3=\\gamma_1=\\gamma_2=1\n",
    "$$\n",
    "\n",
    "データの大きさは`n_sample`とする。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "n_sample = 10_000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# x1, x2, x3, w1, w2, z1, z2の平均\n",
    "rv_mean = [0, 0, 0, 0, 0, 0, 0]\n",
    "\n",
    "# 共分散行列　  x1, x2,x3, w1, w2, z1, z2\n",
    "rv_cov =  [[  1,  0, 0,0.5,  0,  0,  0], # x1\n",
    "           [  0,  1, 0,  0,0.5,  0,  0], # x2\n",
    "           [  0,  0, 1,  0,  0,  0,  0], # x3\n",
    "           [0.5,  0, 0,  1,  0,0.8,  0], # w1\n",
    "           [  0,0.5, 0,  0,  1,  0,0.8], # w2\n",
    "           [  0,  0, 0,0.8,  0,  1,  0], # z1\n",
    "           [  0,  0, 0,  0,0.8,  0,  1]] # z2\n",
    "\n",
    "# ランダム変数の生成\n",
    "rv = multivariate_normal.rvs(rv_mean, rv_cov, size=n_sample)\n",
    "\n",
    "# 説明変数と操作変数の抽出\n",
    "x1 = rv[:,0]  # 外生変数（欠落）\n",
    "x2 = rv[:,1]  # 外生変数（欠落）\n",
    "x3 = rv[:,2]  # 外生変数\n",
    "w1 = rv[:,3]  # 内生変数\n",
    "w2 = rv[:,4]  # 内生変数\n",
    "z1 = rv[:,5]  # 操作変数\n",
    "z2 = rv[:,6]  # 操作変数\n",
    "\n",
    "# yの生成\n",
    "y = 1 + x1 + x2 + x3 +w1 + w2 + np.random.normal(size=n_sample)\n",
    "\n",
    "# DataFrameの作成\n",
    "var_dict = {'y':y,'x1':x1,'x2':x2,'x3':x3,'w1':w1,'w2':w2,'z1':z1,'z2':z2}\n",
    "df_endog2 = pd.DataFrame(var_dict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "上のコードにより，`df_endog2`にシミュレーションで使用するデータが割り当てられた。最初の3行を表示してみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "df_endog2.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 最小二乗法による推定"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "formula_0 = 'y ~ w1 + w2 + x3'\n",
    "print(IV2SLS.from_formula(formula_0, data=df_endog2).fit().summary.tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "欠落変数と相関がない$x_3$の係数の推定値は真の値である`1`に近い。一方，欠落変数と相関する$w_1$と$w_2$は真の値`1`から大きく外れていることが確認できる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 「手計算」"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### 第１段階OLS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "２つの推定式"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "formula_1 = 'w1 ~ z1 + z2 + x3'\n",
    "formula_2 = 'w2 ~ z1 + z2 + x3'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "OLS推定"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "res_endog2_1 = IV2SLS.from_formula(formula_1, data=df_endog2).fit()\n",
    "res_endog2_2 = IV2SLS.from_formula(formula_2, data=df_endog2).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "print(res_endog2_1.summary.tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`z1`の係数の推定値である$\\hat{\\gamma}_1$の$p$値は非常に小さい。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "print(res_endog2_2.summary.tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`z2`の係数の推定値である$\\hat{\\gamma}_2$の$p$値も非常に小さい。\n",
    "\n",
    "次に，それぞれの残差を`df_endog2`に新たな列`w1hat`と`w2hat`として追加する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "df_endog2['w1hat'] = res_endog2_1.fitted_values\n",
    "df_endog2['w2hat'] = res_endog2_2.fitted_values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### 第２段階OLS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "推定式を定義する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "formula_12 = 'y ~ w1hat + w2hat + x3'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "推定結果の表示"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "print(IV2SLS.from_formula(formula_12, data=df_endog2).fit().summary.tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "全ての推定値が真の値である`1`に概ね近いことが分かる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 自動計算"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`linearmodels`の`IV2SLS`を使うが，次の書き方となる。\n",
    "```\n",
    "被説明変数 ~ 定数項 + 外生変数 + [内生変数１+内生変数２ ~ 操作変数１ + 操作変数２]\n",
    "```\n",
    "* 定数項，外生変数がない場合は省いても良い。\n",
    "* 外生変数，操作変数は複数でも可"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "推定式を定義しよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "formula_2 = 'y ~ 1 + x3 + [w1 + w2 ~ z1 + z2]'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "推定結果を表示する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "print(IV2SLS.from_formula(formula_2, data=df_endog2).fit().summary.tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "推定値は「手計算」の結果と同じであることが確認できる。「手計算」の標準誤差とは異なるが，必ず自動計算の値を使うように！"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## シミュレーション：３つの特徴"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 一致性"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "（目的）\n",
    "* `IV`推定量の一致性を確認する。\n",
    "\n",
    "＜シミュレーションの内容＞\n",
    "\n",
    "* 母集団の説明変数には`x1`と`x2`があり，`x2`を欠落変数とする単回帰分析。\n",
    "\n",
    "    $$\n",
    "    y=\\beta_0 + \\beta_1 x_1 + u\n",
    "    $$ (eq:19-8)\n",
    "    \n",
    "* `z`を`x1`操作変数とする。\n",
    "* ２つの推定方法\n",
    "    * 操作変数`z`を使いIV法による推定\n",
    "    * OLSによる推定（この場合，欠落変数バイアスが発生する）\n",
    "* 標本の大きさは$1000$，標本数（ループの回数）を10000として$\\hat{\\beta}_1$の分布を比べる"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "標本数（ループの回数）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "N = 10_000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "母集団のパラメータ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "b0 = 0.5\n",
    "b1 = 1\n",
    "b2 = 0.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "シミュレーションの関数を設定する。\n",
    "* 引数\n",
    "    * `n`：標本の大きさ\n",
    "    * `m`：`x1`と`z`の共分散 $\\text{Cov}(x1,z)=m$\n",
    "    * `ols`（デフォルトは`False`）：OLS推定を一緒に行う場合は`True`\n",
    "* 返り値\n",
    "    * IV推定値のリスト（`ols=False`）\n",
    "    * IV推定値のリストとOLS推定値のリスト（`ols=True`）\n",
    "    \n",
    "（コメント）\n",
    "\n",
    "計算の速度を早めるために下の関数の中では`IV2SLS`は使わず`numpy`の関数を使いIVとOLS推定値を計算している。`IV2SLS`は係数の推定値だけではなく他の多くの統計値も自動的に計算するために一回の計算に比較的に長い時間を要するためである。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def sim_iv(n,m,ols=False):  # n=標本の大きさ, m=x1とｚの共分散\n",
    "    \"\"\"\n",
    "    n: 標本の大きさ\n",
    "    m: 説明変数と対応する操作変数の相関性の度合い\n",
    "    ols: Boolean（デフォルト：False） \n",
    "        False: OLS推定とIV推定も一緒におこなう\n",
    "            返り値：２つのリスト（OLSとIV推定量）\n",
    "        True: IV推定のみ\n",
    "            返り値：１つのリストのみ\n",
    "    \"\"\"\n",
    "    rv_mean = [0, 0, 0]  # x1, x2, z1の平均\n",
    "    rv_cov = [[1.0, 0.3, m],  # x1, x2, z1の共分散行列\n",
    "              [0.3, 1.0, 0.0],  # 全ての変数の分散は１（対角成分）\n",
    "              [m,   0.0, 1.0]]  # Cov(x1,x2)=0.3,  Cov(x2,z)=0, Cov(x1,z)=m,\n",
    "    rv = multivariate_normal.rvs(rv_mean, rv_cov, size=n)  # x1, x2, z1をnセット抽出\n",
    "    x1 = rv[:,0]  # 説明変数\n",
    "    x2 = rv[:,1]  # 欠落変数\n",
    "    z = rv[:,2]  # 操作変数\n",
    "\n",
    "    b1_iv_list = []  # IV推定量を入れる空のリスト\n",
    "    b1_ols_list = []  # OLS推定量を入れる空のリスト\n",
    "    \n",
    "    c = np.ones(n)\n",
    "    \n",
    "    for j in range(N):  # N回のループ\n",
    "        \n",
    "        u = norm.rvs(size=n)  # 母集団の誤差項\n",
    "        y = b0 + b1*x1 + b2*x2 + u  # 母集団回帰式\n",
    "        \n",
    "        # IV 第１ステージ\n",
    "        Xiv1 = np.stack([c,z],axis=1)\n",
    "        pihat = np.linalg.inv((Xiv1.T)@Xiv1)@(Xiv1.T)@x1  # IV推定\n",
    "        x1hat = Xiv1@pihat  # x1の予測値\n",
    "        # IV 第２ステージ\n",
    "        Xiv2 = np.stack([c,x1hat],axis=1)        \n",
    "        beta_iv = np.linalg.inv((Xiv2.T)@Xiv2)@(Xiv2.T)@y  # IV推定\n",
    "        b1_iv_list.append(beta_iv[1])  # IV推定量をリストに追加\n",
    "\n",
    "        if ols==True:  # オプションols=Trueの場合はOLS推定もおこなう\n",
    "            \n",
    "            X = np.stack([c,x1],axis=1)\n",
    "            beta_ols = np.linalg.inv((X.T)@X)@(X.T)@y  # OLS推定\n",
    "            b1_ols_list.append(beta_ols[1])  # OLS推定量をリストに追加\n",
    "\n",
    "        else:  # ols=Falseの場合はOLS推定をおこなわない\n",
    "            pass\n",
    "\n",
    "    if ols==True:  # ols=True の場合の返り値の設定\n",
    "        return b1_iv_list, b1_ols_list\n",
    "    \n",
    "    else:  # ols=False の場合の返り値の設定\n",
    "        return b1_iv_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "シミュレーションの開始"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "b1hat_iv, b1hat_ols = sim_iv(1000,0.9,ols=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "図示"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "xx=np.linspace(0.8,1.3,num=100)  # 図を作成するために横軸の値を設定\n",
    "\n",
    "kde_model_ols=gaussian_kde(b1hat_ols)  # OLS推定量のカーネル密度関数を計算\n",
    "\n",
    "kde_model_iv=gaussian_kde(b1hat_iv)  # IV推定量のカーネル密度関数を計算\n",
    "\n",
    "plt.plot(xx, kde_model_ols(xx), 'g-', label='OLS')  # OLS推定量の分布プロット\n",
    "plt.plot(xx, kde_model_iv(xx),'r-', label='IV')  # IV推定量の分布プロット\n",
    "plt.axvline(x=b1,linestyle='dashed')\n",
    "plt.ylabel('Kernel Density')  # 縦軸のラベル\n",
    "plt.legend()  # 凡例\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "* 欠落変数が発生しているため`OLS`推定量は一致性は満たさない。\n",
    "* `IV`推定量は一致性を満たす。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### 標本の大きさ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "（目的）\n",
    "* `IV`推定量の一致性は大標本特性である。ここでは標本の大きさの効果を確認する。\n",
    "\n",
    "＜シミュレーションの内容＞\n",
    "\n",
    "* 上と同じ設定\n",
    "* `IV`法による単回帰分析のみを考える\n",
    "* 標本の大きさを`30`，`100`，`1000`\n",
    "* それぞれ`N=10000`回推定し$\\hat{\\beta}_1$の分布を比べる"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "上で使った関数`sim_iv()`をデフォルト（`ols=False`）で使う。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "b1hat_iv_30 = sim_iv(30,0.8)\n",
    "b1hat_iv_100 = sim_iv(100,0.8)\n",
    "b1hat_iv_1000 = sim_iv(1000,0.8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "図示"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "xx=np.linspace(0.5,1.5,num=100)  # 図を作成するために横軸の値を設定\n",
    "\n",
    "b1hat_iv_n_list = [b1hat_iv_30, b1hat_iv_100, b1hat_iv_1000]\n",
    "color_list = ['g-', 'r-', 'k-']\n",
    "label_list = ['n=30', 'n=100', 'n=1000']\n",
    "\n",
    "for (i, j, k) in zip(b1hat_iv_n_list, color_list, label_list):\n",
    "    kde_model_iv=gaussian_kde(i)  # IV推定量のカーネル密度推定を設定\n",
    "    b1_dist = kde_model_iv(xx)  # IV推定量のカーネル密度関数を計算\n",
    "    plt.plot(xx, b1_dist, j, label=k)  # IV推定量の分布プロット\n",
    "\n",
    "plt.axvline(x=b1,linestyle='dashed')\n",
    "plt.ylabel('Kernel Density')  # 縦軸のラベル\n",
    "plt.legend()  # 凡例\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "* 標本の大きさが増加すると，分散は低下する。\n",
    "* 標本の大きさが小さい場合，不偏性は満たされないこともわかる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### 操作変数との相関性（弱操作変数）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "（目的）\n",
    "* IV推定法で必須となるのが内生的説明変数と操作変数の相関性である。相関性が高い場合は推定量の標準ごさは低くなるが，逆に相関性が低い場合は推定量の標準誤差が大きくなることを確認する。\n",
    "\n",
    "＜シミュレーションの内容＞\n",
    "\n",
    "* 上と同じ設定\n",
    "* `IV`法による単回帰分析のみを考える\n",
    "* 内生的説明変数と操作変数の共分散`m`を`0.1`，`0.4`，`0.8`の3つのケースを考える。\n",
    "* 標本の大きさ$2000$に固定し，それぞれ`N=10000`回推定し$\\hat{\\beta}_1$の分布を比べる"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "上で使った関数`sim_iv()`を使う。\n",
    "\n",
    "シミュレーションの開始"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "b1hat_iv_weak = sim_iv(2000,0.1)\n",
    "b1hat_iv_mid = sim_iv(2000,0.4)\n",
    "b1hat_iv_strong = sim_iv(2000,0.8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "図示"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "xx=np.linspace(0.5,1.5,num=100)  # 図を作成するために横軸の値を設定\n",
    "\n",
    "b1hat_iv_n_list = [b1hat_iv_weak, b1hat_iv_mid, b1hat_iv_strong]\n",
    "color_list = ['g-', 'r-', 'k-']\n",
    "label_list = ['Cov(x1,z)=0.1', 'Cov(x1,z)=0.4', 'Cov(x1,z)=0.8']\n",
    "\n",
    "for (i, j, k) in zip(b1hat_iv_n_list, color_list, label_list):\n",
    "    kde_model_iv=gaussian_kde(i)  # OLS推定量のカーネル密度推定を設定\n",
    "    b1_dist = kde_model_iv(xx)  # OLS推定量のカーネル密度関数を計算\n",
    "    plt.plot(xx, b1_dist, j, label=k)  # OLS推定量の分布プロット\n",
    "\n",
    "plt.axvline(x=b1,linestyle='dashed')\n",
    "plt.ylabel('Kernel Density')  # 縦軸のラベル\n",
    "plt.legend()  # 凡例\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "$\\text{Cov}(x,z)$が低くなると分散が大きくなり，一致性を失うことも確認できる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## 説明変数の外生性の検定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 説明"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "次式を考えよう。\n",
    "\n",
    "$$\n",
    "y=\\beta_0+\\beta_1w+u\n",
    "$$ (eq:19-9)\n",
    "\n",
    "* $\\text{Cov}(w,u)=0$の場合，`OLS`がより良い推定量\n",
    "    * `OLS`推定量は不偏性・一致性を満たす。\n",
    "    * `IV-2SLS`推定量一致性を満たすが不偏性は欠き，標準誤差が大きくなり正確性を損なう。\n",
    "* $\\text{Cov}(w,u)\\neq 0$の場合，`IV-2SLS`がより良い推定量\n",
    "    * `OLS`推定量は不偏性・一致性を満たさない。\n",
    "    * `IV-2SLS`推定量は一致性を満たす。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "このように説明変数が外生的か内生的かによって推定量の性質が大きく異なる。では`IV`推定法と`OLS`推定法のどちらを使うべきなのか。この問題は説明変数が外生的か（$\\text{Cov}(w,u)=0$）それとも内生的か（$\\text{Cov}(w,u)\\neq 0$）に依存しており，ここではその検定方法について説明する。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "上で考えたケース２を使って確かめる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 検定方法１：`Wu-Hausman`検定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "帰無仮説と対立仮説：\n",
    "\n",
    "$\\text{H}_0:\\;w$は外生変数（$\\text{Cov}(w,u)=0$）\n",
    "\n",
    "$\\text{H}_0:\\;w$は内生変数（$\\text{Cov}(w,u)\\neq0$）\n",
    "\n",
    "\n",
    "---\n",
    "（基本的な考え方）\n",
    "\n",
    "帰無仮説のもとでは，`IV`推定法も`OLS`推定法も一致性も満たすため\n",
    "\n",
    "$$\\left|\\hat{\\beta}_j^{\\text{OLS}}-\\hat{\\beta}_j^{\\text{IV}}\\right|$$\n",
    "\n",
    "の値は小さいはず。逆に，帰無仮説が成立しない場合，`OLS`推定量は不偏性も一致性も失うことになり，上の値は大きくなる。即ち，上の値が大きければ（小さければ），帰無仮説を棄却できる可能性が高くなる（低くなる）。この考えを利用したのが`Wu-Hausman`検定である。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "`res_2`のメソッド`.wu_hausman()`を使うと検定統計量と$p$値が表示される。\n",
    "\n",
    "（注意）`.wu_hausman()`であって`.wu_hausman`ではない。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "res_2.wu_hausman()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "$p$値は0.0948。5%の有意水準では帰無仮説を棄却できない（外生性を棄却できない）が，10%では棄却できる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 検定方法２"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "まず検定方法について説明し，「考え方」については後述する。次の回帰式を考える。\n",
    "\n",
    "$$\n",
    "y=\\beta_0+\\beta_1x + \\beta_2w+u\n",
    "$$ (eq:19-10)\n",
    "\n",
    "* $x$は外生変数\n",
    "* $w$は内生性が疑われる説明変数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "２段階で検定する。\n",
    "\n",
    "**第１段階**\n",
    "* 式[](eq:19-10)のOLS残差$\\hat{u}$を計算する。\n",
    "\n",
    "**第２段階**\n",
    "* 式[](eq:19-10)に$\\hat{u}$を加えてOLS推定する。\n",
    "\n",
    "    $$y=\\gamma_0+\\gamma_1x + \\gamma_2w+\\gamma_u\\hat{u}+e\\qquad\\qquad (\\text{式２})$$\n",
    "\n",
    "* 次の検定をおこなう。\n",
    "    * $\\text{H}_0:\\;\\hat{\\gamma}_u=0$（$w$は外生変数である）\n",
    "    * $\\text{H}_A:\\;\\hat{\\gamma}_u\\neq 0$（$w$は内生変数である）\n",
    "\n",
    "（コメント）計算上$\\hat{\\gamma}_i=\\hat{\\beta}_i,\\;i=0,1,2$が成り立つことになる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "**-------- 考え方 --------**\n",
    "\n",
    "＜第１段階＞\n",
    "* (式１)を推定すると$y$を$\\hat{y}$と$\\hat{u}$に分解することができる。\n",
    "\n",
    "    $$y=\\hat{y}+\\hat{u}$$\n",
    "    \n",
    "    * $\\hat{y}$：説明変数で説明できる$y$の部分\n",
    "    * $\\hat{u}$：$y$のその他の部分\n",
    "* $u$が$w$と相関していれば，その相関する部分が$\\hat{u}$に含まれることになる。以下のようなイメージ。\n",
    "\n",
    "    $$\\hat{u}=f\\left(\\hat{u}_w,v\\right)$$\n",
    "    \n",
    "    * $\\hat{u}_w$：$w$と相関する部分\n",
    "    * $v$：$w$と相関しない部分\n",
    "\n",
    "＜第２段階＞\n",
    "* (式２)を推定すると，相関部分である$\\hat{u}_w$を$\\gamma_u\\hat{u}$として取り出すことができる。\n",
    "* もし$\\gamma_u$が有意であれば，$\\text{H}_0$は棄却される。$\\hat{u}_w$は存在するということになり，$w$は内生変数と判断できる。\n",
    "* もし$\\gamma_u$が有意でなければ，$\\text{H}_0$は棄却できない。$\\hat{u}_w$は存在しないということになり，$w$は外生変数と判断できる。"
   ]
  },
  {
   "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": [
    "form_2a = 'educ ~ 1 + exper + expersq + motheduc + fatheduc'\n",
    "\n",
    "mod_2a = IV2SLS.from_formula(form_2a, data=mroz)\n",
    "\n",
    "res_2a = mod_2a.fit(cov_type='unadjusted')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "**第２段階**\n",
    "\n",
    "`res_2a`の残差を使うが，メソッド`resids`を使い取得する。すなわち，`res_2a.resids`が残差を返す。また次の点を覚えておこう。\n",
    "* `res_2a.resids`を直接回帰式に入れることはできないが，`I()`の引数とすれば可能となる。これは`statsmodels`と同じである。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "form_2b = 'lwage ~ 1 + educ + exper + expersq + I(res_2a.resids)'\n",
    "\n",
    "mod_2b = IV2SLS.from_formula(form_2b, data=mroz)\n",
    "\n",
    "res_2b = mod_2b.fit(cov_type='unadjusted')\n",
    "\n",
    "print(res_2b.summary.tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "$p$値は0.0928で，5%の有意水準では帰無仮説を棄却できない（内生性を棄却できない）。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### 自動計算"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`linearmodels`には上の計算を自動でおこなった結果を示すメソッドが用意されている。`res_2`の`.wooldridge_regression`である。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "res_2.wooldridge_regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "（注意）「手計算」で使った`.fit()`のオプション`debiased=True`を使うと異なる数値になる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## 操作変数の妥当性検定（過剰識別制約検定）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 説明"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "次の回帰式を考えよう。\n",
    "\n",
    "$$\n",
    "y=\\beta_0+\\beta_1x+\\beta_2w+u\n",
    "$$ (eq:19-11)\n",
    "\n",
    "* $x$：外生変数\n",
    "* $w$：内生性が疑われる説明変数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "**＜ケース１：１つの内生的説明変数に１つの操作変数がある場合＞**\n",
    "* 操作変数$z$\n",
    "\n",
    "この場合，操作変数の外生性$\\text{Cov}(z,u)=0$が満たされないといけないが，これを検定できない。\n",
    "\n",
    "（理由）\n",
    "* 式[](eq:19-11)の`w`は内生性が疑われるため，その式を使い計算したOLS残差$\\hat{u}$は真の誤差項を捉えていない可能性がある。従って，この$\\hat{u}$を使っても意味がある検定とはならない。\n",
    "* 式[](eq:19-11)に$z$を使い残差を計算することも考えられるが，そもそも$z$の妥当性が分からないため，これも真の誤差項を捉えていない可能性がある。\n",
    "\n",
    "この場合，経済理論に基づいて操作変数の外生性を正当化できるかが問題になる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "**＜ケース２：１つの内生的説明変数に複数の操作変数がある場合＞**\n",
    "* 例えば，２つの操作変数$z_1$と$z_2$\n",
    "\n",
    "この場合，**同時**に$\\text{Cov}(z_1,u)=\\text{Cov}(z_2,u)=0$が成立するかを検定する方法があり，その１つがSargan検定と呼ばれる。\n",
    "\n",
    "（注意）\n",
    "* Sargan検定は，全ての操作変数（上の例では２つ）が**同時**に妥当かどうかを検定する。言い換えると，同時に妥当ではない場合は，次のパターンとなる。\n",
    "    1. $z_1$は妥当だが，$z_2$は妥当ではない。\n",
    "    1. $z_1$は妥当ではないが，$z_2$は妥当である。\n",
    "    1. $z_1$と$z_2$の両変数は妥当ではない。\n",
    "    \n",
    "従って，同時に妥当ではない場合，どの操作変数が有効かを調べることはできない。\n",
    "\n",
    "（コメント）\n",
    "* 操作変数の数が内生変数の数を上回っているケースは「過剰識別（overidentified）」と呼ばれる。\n",
    "* 操作変数の数と内生変数の数が等しい場合は「適度識別」（just identified）と呼ばれる。\n",
    "\n",
    "Sargan検定は過剰識別の場合のみ意味がある検定であり，「過剰識別制約検定」の一つとなる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "**Sargan検定**\n",
    "\n",
    "２段階で検定する。\n",
    "\n",
    "**第１段階**\n",
    "* IV/2SLS推定法で推定し，残差$\\hat{u}$を取得する。\n",
    "\n",
    "**第２段階**\n",
    "* 帰無仮説と対立仮説の設定\n",
    "    * $\\text{H}_0:$ 全ての操作変数は妥当である（$\\hat{u}$は全ての操作変数と相関性がなく外生的である）。\n",
    "    * $\\text{H}_A:$ $\\text{H}_0$は成立しない（少なくとも１つの操作変数は$\\hat{u}$と相関する）。\n",
    "    * 上の例では，操作変数の数は`2`であり内生変数の数は`1`であるため$q=2-1=1$の過剰識別制約があると考える。\n",
    "* $\\hat{u}$を被説明変数として全ての外生変数と操作変数に対してOLS推定し，決定係数$R^2$を計算する。\n",
    "* $LM$統計量 $=nR^2$\n",
    "    * $nR^2\\sim\\chi(q)$\n",
    "    * $q>1$は操作変数の数と内生的説明変数の差\n",
    "    * $n$は標本の大きさ\n",
    "\n",
    "（帰無仮説棄却の意味）\n",
    "\n",
    "* 少なくとも１つの操作変数は誤差項と相関性があるということになる。（欠落変数とも解釈可能）\n",
    "* しかしどの操作変数が内生的かは分からない。\n",
    "\n",
    "---\n",
    "上で考えたケース２を使って確かめる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 「手計算」"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "第１段階の計算結果として`res_2`を使う。\n",
    "\n",
    "第２段階の計算のために`res_2`の属性`.resids`を使い回帰残差を取得し，直接以下の回帰式に入れる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "form_aux = 'I(res_2.resids) ~ 1 + motheduc + fatheduc + exper + expersq' # 外生的説明変数を省いてもよい\n",
    "\n",
    "mod_aux = IV2SLS.from_formula(form_aux, data=mroz)\n",
    "\n",
    "res_aux = mod_aux.fit(cov_type='unadjusted')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "検定統計量の計算"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "r2 = res_aux.rsquared  # 決定係数\n",
    "n = res_aux.nobs    # 標本の大きさ\n",
    "teststat = n * r2      # 検定統計量\n",
    "pval = 1 - chi2.cdf(teststat, 1)  # p値の計算\n",
    "print(pval)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "５％有意水準では帰無仮説を棄却できない。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 自動計算 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`res_2`のメソッド`.sargan`を使うと上と同じ計算結果を表示できる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "res_2.sargan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## 同時方程式モデルと`IV`推定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### 同時性バイアス"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "同時方程式モデルとは，均衡メカニズムなどを通して複数の内生変数が複数の式によって同時決定されるモデルである。例として労働の需給モデルを考えよう。均衡では需要量（$L_d$）と供給量（$L_s$）は等しくなり（$L=L_d=L_s$），需要と供給はそれぞれ均衡賃金（$W$）に依存する。\n",
    "* 供給関数\n",
    "\n",
    "    $$\n",
    "    L = s_0+s_1 W + s_2 X_s + u_s\n",
    "    $$ (eq:19-12)\n",
    "\n",
    "    * $s_1>0$\n",
    "    * $X_s=$ 供給の「その他」の決定要因（例えば，限界費用）\n",
    "    * $u_s=$ 供給の誤差項\n",
    "\n",
    "* 需要関数\n",
    "\n",
    "    $$\n",
    "    W = d_0+d_1 L + d_2 X_d + u_d\n",
    "    $$ (eq:19-13)\n",
    "    \n",
    "    * $d_1<0$\n",
    "    * $X_d=$ 需要の「その他」の決定要因（例えば，所得）\n",
    "    * $u_d=$ 需要の誤差項\n",
    "\n",
    "（相関性の仮定）\n",
    "* $\\text{Cov}(X_s,u_s)=\\text{Cov}(X_s,u_d)=0$\n",
    "* $\\text{Cov}(X_d,u_d)=\\text{Cov}(X_d,u_s)=0$\n",
    "* $\\text{Cov}(u_s,u_d)=0$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### `IV`推定法の適用"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "同時性バイアスは`IV`推定法で対処可能である。考え方は簡単である。\n",
    "\n",
    "供給曲線式[](eq:19-12)の推定\n",
    "* $W$の操作変数として$X_d$を使う。$X_d$は操作変数の３つの条件を満たす。\n",
    "    * 式[](eq:19-12)にない\n",
    "    * $\\text{Cov}(X_d,u_s)=0$（相関性の仮定から）\n",
    "    * $\\text{Cov}(X_d,W)\\neq 0$（式[](eq:19-13)から）\n",
    "\n",
    "供給曲線式[](eq:19-13)の推定\n",
    "* $L$の操作変数として$X_s$を使う。$X_s$は操作変数の３つの条件を満たす。\n",
    "    * 式[](eq:19-13)にない\n",
    "    * $\\text{Cov}(X_s,u_d)=0$（相関性の仮定から）\n",
    "    * $\\text{Cov}(X_s,L)\\neq 0$（式[](eq:19-12)から）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 例"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### データ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "データセット`mroz`を使う。\n",
    "* 労働供給曲線：労働市場に参加する既婚女性の労働供給関数\n",
    "* 労働需要曲線：企業が提示（オファー）する賃金をその決定要因の関数として表す\n",
    "\n",
    "上の記号に対応する変数をリストアップする（1975年のデータ）。\n",
    "* $L$：`hours`（労働時間）\n",
    "* $W$：`lwage`（賃金時間額の対数）\n",
    "* $X_s$：労働供給の外生変数\n",
    "    * `age`（年齢）\n",
    "    * `kidslt6`（子どもが６歳未満）\n",
    "    * `nwifeinc`：（(家計の所得 - 賃金\\*時間)/1000）\n",
    "* $X_d$：賃金オファーの外生変数\n",
    "    * `exper`（雇用経験）\n",
    "    * `expersq`（雇用経験の２乗）\n",
    "* 両方に含まれる変数：`educ`（教育年数）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### 供給曲線の推定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "**`IV`推定**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "form_L = 'hours ~ 1 + educ + age + kidslt6 + nwifeinc + [lwage ~ exper + expersq]'\n",
    "\n",
    "mod_L = IV2SLS.from_formula(form_L, data=mroz)\n",
    "\n",
    "res_L = mod_L.fit(cov_type='unadjusted')\n",
    "\n",
    "print(res_L.summary.tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "**`OLS`推定**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "form_L_ols = 'hours ~ 1 + educ + age + kidslt6 + nwifeinc + lwage '\n",
    "\n",
    "mod_L_ols = IV2SLS.from_formula(form_L_ols, data=mroz)\n",
    "\n",
    "res_L_ols = mod_L_ols.fit(cov_type='unadjusted')\n",
    "\n",
    "print(res_L_ols.summary.tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "２つの推定結果を比べると，推定値を含めて大きく異なることがわかる。２段回推定法により同時性バイアスを取り除いた結果と解釈できる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### 需要曲線の推定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "**`IV`推定**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "form_P = 'lwage ~ 1 + educ + exper + expersq + [hours ~ age + kidslt6 + nwifeinc]'\n",
    "\n",
    "mod_P = IV2SLS.from_formula(form_P, data=mroz)\n",
    "\n",
    "res_P = mod_P.fit(cov_type='unadjusted')\n",
    "\n",
    "print(res_P.summary.tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "**`OLS`推定**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "form_W_ols = 'lwage ~ 1 + educ + exper + expersq + hours'\n",
    "\n",
    "mod_W_ols = IV2SLS.from_formula(form_W_ols, data=mroz)\n",
    "\n",
    "res_W_ols = mod_W_ols.fit(cov_type='unadjusted')\n",
    "\n",
    "print(res_W_ols.summary.tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "労働供給の推定結果と比べて，パラメータの推定値の変化は大きくないが，統計的優位性は大きく変化している。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## 測定誤差と`IV`推定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### 測定誤差によるバイアス"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "ここではシミュレーションを通して，測定誤差バイアスがある場合でもIV推定法を使うことで一致性を満たす推定量を得ることが可能であることを示す。\n",
    "\n",
    "次の母集団回帰式を考えよう。\n",
    "\n",
    "$$\n",
    "y=\\beta_0+\\beta_1 x^* + \\eta\n",
    "$$ (eq:19-14)\n",
    "\n",
    "ここで，被説明変数$y$は正確に測定できるが，説明変数$x^*$は以下の式に従って測定される仮定する。\n",
    "\n",
    "$$\n",
    "x=x^*+e\n",
    "$$ (eq:19-15)\n",
    "\n",
    "* $x$：測定値\n",
    "* $e$：測定誤差\n",
    "\n",
    "（仮定）\n",
    "\n",
    "* 測定誤差$e$は真の値と無関係。即ち，$\\text{Cov}(x^*,e)=0$\n",
    "\n",
    "（結果）\n",
    "\n",
    "* 次式をOLS推定する場合，$\\hat{\\beta}_1$は不偏性・一致性を満たさない。\n",
    "\n",
    "$$\n",
    "y=\\beta_0+\\beta_1 x + u,\\qquad u=\\eta-\\beta e\n",
    "$$ (eq:19-16)\n",
    "\n",
    "（理由）\n",
    "\n",
    "仮定４：$\\text{Cov}(x,u)=0$が満たされない。\n",
    "\n",
    "$$\n",
    "\\text{Cov}(x,u)=\\text{E}[xu]\n",
    "=\\text{E}\\left[(x^*+e)(\\eta-\\beta e)\\right]\n",
    "=-\\beta\\cdot\\text{E}(e^2)>0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### `IV`推定法の適用"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`IV`推定法での対処方法を考える。$x^*$を測定する変数$z$があるとする。\n",
    "\n",
    "$$\n",
    "z=x^*+v\n",
    "$$ (eq:19-17)\n",
    "\n",
    "$v\\sim\\text{iid}(0,\\sigma)$は誤差であり，$\\text{iid}$は独立同一分布を意味する。\n",
    "\n",
    "**`IV`推定**\n",
    "* $x$の操作変数として$z$を使う\n",
    "* 以下の条件を満たすと仮定する。\n",
    "    1. もとの式に含まれていない。\n",
    "    1. $\\text{Cov}(v,u)=0\\;\\;\\Rightarrow\\;\\;\\text{Cov}(z,u)=0$\n",
    "    1. $\\text{Cov}(z,x)\\neq 0$（$z$も$x$も$x^*$を測定する変数）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### シミュレーション"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "（目的）\n",
    "* 測定誤差によるバイアスを示す。\n",
    "* `IV`推定法により一致性が成立することを示す。\n",
    "\n",
    "＜シミュレーションの内容＞\n",
    "\n",
    "* 単回帰分析\n",
    "\n",
    "    $$\n",
    "    y=\\beta_0 + \\beta_1 x + u\n",
    "    $$ (eq:19-18)\n",
    "    \n",
    "* ２つのケース\n",
    "    1. $OLS$推定\n",
    "    2. $IV$推定\n",
    "* それぞれのケースで標本の大きさ$n=100$\n",
    "* 1000回推定し$\\hat{\\beta}_1$の分布を比べる"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "標本の大きさと標本数（ループの回数）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "n = 1000\n",
    "N = 10_000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "母集団のパラメータの真の値"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "b0 = 1.0\n",
    "b1 = 1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "x_pop = uniform.rvs(1,10,size=n)  # 母集団の説明変数\n",
    "u = norm.rvs(scale=1, size=n)  # 母集団の誤差項\n",
    "y = b0 + b1*x_pop + u  # 母集団回帰式"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "測定誤差の標準偏差"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "error_sd = 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "シミュレーション開始"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# シミュレーションで計算した推定量を入れる空のリストの作成\n",
    "b1_ols_list = []  # OLS推定量\n",
    "b1_iv_list = []  # IV推定量\n",
    "\n",
    "for j in range(N):  # N回のループ\n",
    "    x = x_pop + norm.rvs(scale=error_sd, size=n)  # 測定誤差\n",
    "    z = x_pop + norm.rvs(scale=error_sd, size=n)  # 操作変数\n",
    "    \n",
    "    c = np.ones(n)  # 定数項\n",
    "    \n",
    "    # IV 第１ステージ\n",
    "    Xiv1 = np.stack([c,z],axis=1)\n",
    "    pihat = np.linalg.inv((Xiv1.T)@Xiv1)@(Xiv1.T)@x  # IV推定\n",
    "    xhat = Xiv1@pihat  # x1の予測値\n",
    "    # IV 第２ステージ\n",
    "    Xiv2 = np.stack([c,xhat],axis=1)        \n",
    "    beta_iv = np.linalg.inv((Xiv2.T)@Xiv2)@(Xiv2.T)@y  # IV推定\n",
    "    b1_iv_list.append(beta_iv[1])  # b1のIV推定量をリストに追加\n",
    "\n",
    "    # OLS\n",
    "    X = np.stack([c,x],axis=1)        \n",
    "    beta_ols = np.linalg.inv((X.T)@X)@(X.T)@y  # OLS推定\n",
    "    b1_ols_list.append(beta_ols[1])  # b1のOLS推定量"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "結果の図示"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "xx=np.linspace(0.4,1.2,num=100)  # 図を作成するために横軸の値を設定\n",
    "\n",
    "kde_model_ols=gaussian_kde(b1_ols_list)  # t値のカーネル密度推定を設定\n",
    "b1_ols_dist = kde_model_ols(xx)\n",
    "\n",
    "kde_model_iv=gaussian_kde(b1_iv_list)  # t値のカーネル密度推定を設定\n",
    "b1_iv_dist = kde_model_iv(xx)\n",
    "\n",
    "plt.plot(xx, b1_ols_dist, 'g-', label='OLS Estimates')  # t値の分布プロット\n",
    "plt.plot(xx, b1_iv_dist,'r-', label='IV Estimates')  # t分布\n",
    "plt.axvline(x=b1,linestyle='dashed')\n",
    "plt.ylabel('Kernel Density')  # 縦軸のラベル\n",
    "plt.legend()\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "* OLS推定量は不偏性も一致性も満たさない。\n",
    "* IV推定量は一致性を満たす。"
   ]
  }
 ],
 "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.5"
  },
  "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": 2
}