21_TruncregTobitHeckit.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 japanize_matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.formula.api as smf\n",
    "import py4macro\n",
    "import wooldridge\n",
    "\n",
    "from py4etrics.truncreg import Truncreg\n",
    "from py4etrics.tobit import Tobit\n",
    "from py4etrics.heckit import Heckit\n",
    "from py4etrics.hetero_test import het_test_probit\n",
    "from scipy.stats import norm\n",
    "\n",
    "# 警告メッセージを非表示\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## 説明"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "**制限従属変数**（Limited Dependent Variables; LDV）とは，被説明変数が取り得る値が大きく制限される場合を指す。`Logit`と`Probit`も制限従属変数の例であり，被説明変数は$(0,1)$に制限されている。この章では，違ったタイプの制限従属変数を考える。\n",
    "\n",
    "データの性質上２つに分けることができる。\n",
    "1. 切断データ（Truncated Data）\n",
    "    * 特定のデータが標本に含まれていない。\n",
    "    * 例：平成30年度「環境にやさしい企業行動調査」の調査対象企業は，従業員数500名以上となっており，500名未満の企業は含まれていない。\n",
    "1. 打ち切りデータ（Censored Data）\n",
    "    1. 端点解の場合\n",
    "        * 募金の金額の例：募金しない場合は０円だが，募金する人の額は正の値を取る。（下からの打ち切りデータ）\n",
    "    1. データ制約の場合\n",
    "        * 所得調査の例：Ｘ万円からＹ万円のようにカテゴリーでまとめる場合が通常であり，最高額カテゴリーはＺ万円以上のようになる。この場合，Ｚ万円以上の所得は全てＺ万円となる。（上からの打ち切りデータ）\n",
    "\n",
    "（コメント）以下のようにも呼ばれる\n",
    "* **下から**の打ち切りデータ ＝ **左から**の打ち切りデータ（left-censored）\n",
    "* **上から**の打ち切りデータ ＝ **右から**の打ち切りデータ（right-censored）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "データの性質に分けて次のモデルの使い方を説明する。\n",
    "* 切断回帰モデル（Truncated Regression）：切断データ\n",
    "* Tobitモデル：打ち切りデータ\n",
    "* Heckitモデル：切断データで選択バイアス問題が発生する場合"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## 切断回帰モデル"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### 説明"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "無作為な形ではなく，ある特定の一部のデータが標本から欠落している切断データの分析を考える。例として，女性の労働供給を考えよう。\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "&y_i=\\beta_0+\\beta_1x_i+u_i\\qquad i\\in\\cal{Y}\\\\\n",
    "&(y_i,x_i)\\text{は観察されない}\\qquad i\\in\\cal{N}\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "* $\\cal{Y}$：観測されたデータの集合（通常，無作為に選択されないと考える）。\n",
    "* $\\cal{N}$：観測されないデータの集合（通常，無作為に選択されないと考える）。\n",
    "* $x$：労働供給に関する決定要因（例えば，教育年数）\n",
    "* $u|x\\sim \\text{Normal}(0,\\sigma^2)$：労働供給に関するランダムな要素（例えば，好み）\n",
    "    * $x$を所与とすると，$u$は正規分布に従う。\n",
    "    * この仮定が非常に重要であり計算の前提となる。\n",
    "* $y_i$：実際の労働供給（非負の連続変数）\n",
    "\n",
    "＜結果＞\n",
    "* 標本集合$\\cal{Y}$は無作為標本ではない。従って，GM仮定２が満たされないためOLS推定量$\\hat{\\beta}_{\\text{ols}}$は不偏性を満たさない。また一致性も満たさない。\n",
    "\n",
    "---\n",
    "まずシミュレーションを使ってこの結果を直感的に確認する。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### シミュレーション"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# 標本の大きさ \n",
    "n = 100\n",
    "\n",
    "# y*を決定するx\n",
    "x = np.sort(norm.rvs(0,3,size=n))  # ランダム変数を生成し昇順に並べる\n",
    "\n",
    "# 被説明変数\n",
    "y = 1 + x + norm.rvs(0,3,size=n)\n",
    "\n",
    "# DataFrame\n",
    "df = pd.DataFrame({'Y':y, 'Y_trunc':y, 'X':x})\n",
    "\n",
    "# 閾値\n",
    "left = 0\n",
    "\n",
    "# 切断データの作成\n",
    "cond = (df.loc[:,'Y'] <= left)\n",
    "df.loc[cond,'Y_trunc'] = np.nan\n",
    "\n",
    "# 切断された被説明変数\n",
    "y_trunc = df['Y_trunc']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "ここで`np.nan`とは`NumPy`の定数であり，欠損値`NaN`を示す。`cond`の条件にある閾値`left=0`以下のデータは欠損値として設定されている。`NumPy`には次の定数もあるのでこの機会に紹介する。\n",
    "* `np.nan`：`NaN`\n",
    "* `np.inf`：（正の）無限\n",
    "* `np.e`：$e=2.71...$\n",
    "* `np.pi`：$\\pi=3.14...$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "母集団回帰式"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "formula_full = 'Y ~ X'\n",
    "\n",
    "result_full=  smf.ols(formula_full, data=df).fit()\n",
    "\n",
    "b0_full,b1_full = result_full.params\n",
    "\n",
    "result_full.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "切断データを使ったOLS回帰"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "formula_trunc = 'Y_trunc ~ X'\n",
    "\n",
    "result_trunc = smf.ols(formula_trunc, data=df).fit()\n",
    "\n",
    "b0_trunc, b1_trunc = result_trunc.params\n",
    "\n",
    "result_trunc.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "比較するために２つの結果を図示する。同じ図に描くために，先に切断データを整理する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "x_trunc = df.dropna(subset=['Y_trunc']).loc[:,'X']  #1\n",
    "x_min = min(x_trunc)  #2\n",
    "x_max = max(x_trunc)  #3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "* `#1`：`.dropna()`を使って列`Y_trunc`から欠損値がある行を取り除き，列`X`を`x_trun`に割り当てる。\n",
    "* `#2`：`x_trun`の最小値を`x_min`に割り当てる。\n",
    "* `#3`：`x_trun`の最大値を`x_max`に割り当てる。\n",
    "\n",
    "`matplotlib`を使って２つの図を重ねて図示してみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "#1：母集団データの散布図\n",
    "plt.scatter(x, y, facecolors='none', edgecolors='gray',label=r'$y^{*}$')\n",
    "\n",
    "#2：切断データの散布図\n",
    "plt.scatter(x, y_trunc, facecolors='gray', label=r'$y$ and $y^{*}$')\n",
    "\n",
    "# 母集団OLS\n",
    "plt.plot(x, b0_full+b1_full*x, 'k', linewidth=3, label='Population: $y^*$')\n",
    "\n",
    "# 切断回帰\n",
    "plt.plot(x_trunc, b0_trunc+b1_trunc*x_trunc, 'r', lw=3,label=r'Sample: $y> 0$')\n",
    "\n",
    "plt.xlabel('x')\n",
    "plt.ylabel(r'$y$ and $y^{*}$')\n",
    "plt.legend()\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "＜`matplotlib`コードの説明＞\n",
    "\n",
    "* `#1`：母集団データの散布図\n",
    "    * 引数`facecolors`は全ての観測値の点の色を「なし」に設定する。\n",
    "    * 引数`edgecolors`は全ての観測値の円周をグレーに設定する。\n",
    "    * この２つにより，全ての観測値は左下の観測値のように白抜きのグレーの円として表示される。\n",
    "* `#1`：切断データの散布図\n",
    "    * 後にくる`plt.plot()`は前にある`plt.plot()`を上書きする。従って，切断データの観測値は母集団データの散布図上に重ねて表示されるため，引数`facecolors`を`gray`に設定することにより切断データだけをグレーの点として表示させることができる。"
   ]
  },
  {
   "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_1 x+e\\qquad\\text{(式１)}$$\n",
    "\n",
    "であり，全ての観測値を使う（即ち，切断されたデータがない）場合の条件付き期待値である母集団回帰関数（母集団で平均で成立する回帰線）は、$\\text{E}(e)=0$の仮定の下では，次式で与えられる。\n",
    "\n",
    "$$\\text{E}\\left(y|x\\right)=\\beta_0+\\beta_1 x\\qquad\\text{(式２)}$$\n",
    "\n",
    "逆に言うと，(式２)に平均0の誤差項$e$を加えたものが母集団回帰式である。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "**＜切断回帰＞**\n",
    "\n",
    "$y>0$のデータだけを使う回帰式は\n",
    "\n",
    "$$y_i=\\beta_0+\\beta_1 x_i+u_i\\qquad\\text{(式３)}$$\n",
    "\n",
    "である。しかし，標本では$y_i\\leq 0$のデータが切断されているため，$y_i>0$を考慮し$x_i$を所与として(式３)の期待値を取ると\n",
    "\n",
    "$$\\text{E}\\left(y_i|y_i>0,x_i\\right)=\\beta_0+\\beta_1 x_i+\\sigma_u\\lambda\\left(w_i\\right)\\qquad\\text{(式４)}$$\n",
    "\n",
    "となる。ここで\n",
    "* $w_i\\equiv\\dfrac{A-(\\beta_0+\\beta_1x_i)}{\\sigma}$\n",
    "* $A$：下限の値（この例では`0`）\n",
    "* $\\sigma_u$：誤差項の標準偏差\n",
    "\n",
    "である。$\\sigma_u\\lambda\\left(w_i\\right)$の導出や$\\lambda\\left(\\cdot\\right)$の関数形が重要ではない。重要なのは，逆ミルズ比（inverse Mill's ratio）と呼ばれる$\\lambda(w_i)\\neq 0$の存在であり，$\\lambda(w_i)$は$y>0$となる確率の影響を捉えているという点である。従って，$\\beta_0$と$\\beta_1$の推定には，(式４)に残差$\\nu_i$を追加した次式\n",
    "\n",
    "$$\n",
    "y_i=\\beta_0+\\beta_1 x_i+\\sigma_u\\lambda\\left(w_i\\right)+\\nu_i\\qquad\\text{(式５)}\n",
    "$$\n",
    "\n",
    "を推定する必要がある。しかし(式３)をOLS推定すると，$\\lambda(x)$が欠落することになり，**欠落変数バイアス**が発生することになる。このバイアスは，上の図で黒と赤の回帰線の傾きの違いに現れている。OLS推定量はゼロ方向にバイアスが発生する事が知られている。この問題を克服するのが切断回帰モデルである。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### `Truncreg`モジュールの使い方"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "切断データを扱うために切断回帰モデルを展開する。ここでは具体的な内容は割愛するが，`Logit`と`Probit`と同じように最尤法をつかい推定する。\n",
    "\n",
    "* 切断回帰モデルの推定量は一致性を満たす。\n",
    "\n",
    "---\n",
    "`statsmodels`も`linearmodels`も`Tobit`推定のモジュールがない。その代わりに著者が作成した`py4etrics`パッケージの関数`trucreg`モジュールを使い推定する。これは`statsmodels`の`GenericMaximumLikelihoodModel`を使い実装したものである。使用する上で[このサイト](https://www.statsmodels.org/devel/examples/notebooks/generated/generic_mle.html)にある次の点に注意する必要がある。\n",
    "* `R`の推定値と小数点第４位まで同じになるが，標準偏差は小数点第２位までが同じとなる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "````{note}\n",
    "MacではTerminal、WindowsではGit Bashを使い、次のコマンドで`py4etrics`モジュールをインストールできる。\n",
    "```\n",
    "pip install git+https://github.com/spring-haru/py4etrics.git\n",
    "```\n",
    "````"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true,
    "tags": [
     "output_scroll"
    ]
   },
   "source": [
    "\n",
    "---\n",
    "＜使い方＞\n",
    "\n",
    "基本的に`statsmodels`の`ols`と同じだが，追加的な操作とオプションがある。\n",
    "1. 推定式を決める\n",
    "```\n",
    "    formula = 'y ~ 1 + x'\n",
    "```\n",
    "2. `Truncreg`の`from_formula`モジュールを使って推定\n",
    "```\n",
    "    Truncreg.from_formula(formula, left=<A>, right=<B>, data=<C>).fit()\n",
    "```\n",
    "\n",
    "ここで \n",
    "* `left`:左切断の値（デフォルトは$-\\infty$であり，「切断なし」という意味）\n",
    "* `right`:右切断の値（デフォルトは$\\infty$，「切断なし」という意味）\n",
    "* `deta`:データの指定\n",
    "* 切断方向の設定：\n",
    "    * `left`だけに値を設定する場合は左切断回帰（left-truncated）となる。\n",
    "    * `right`だけに値を設定する場合は右切断回帰（right-truncated）となる。\n",
    "    * `left`と`right`の両方に値を設定する場合は左右切断回帰（left- and right-truncated）となる。\n",
    "    * `left`と`right`の両方に値を設定しない場合は通常の最尤推定となる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### 切断回帰推定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "例として`wooldridge`の`HTV`のデータを使い推定する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true,
    "tags": [
     "output_scroll"
    ]
   },
   "outputs": [],
   "source": [
    "htv = wooldridge.data('htv')\n",
    "wooldridge.data('htv', description=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "＜目的＞\n",
    "\n",
    "教育（`educ`）が賃金（`wage`）に与える影響を探る。1991年の時間賃金を対数化した`lwage`を被説明変数として使い，次の説明変数を使う。\n",
    "* `educ`：1991年までに修了した最高学位の指標\n",
    "* `abil`：能力を捉える指標\n",
    "* `exper`：潜在的な労働経験\n",
    "* `nc`：米国北中部のダミー変数\n",
    "* `west`：米国西部のダミー変数\n",
    "* `south`：米国南部のダミー変数\n",
    "* `urban`：都市部のダミー変数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "まずOLS推定を行う。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "formula_trunc = 'lwage ~ 1 + educ + abil + exper + nc + west + south + urban'\n",
    "\n",
    "res_ols = smf.ols(formula_trunc, data=htv).fit()\n",
    "\n",
    "print(res_ols.summary().tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`educ`の係数は`0.1037`であり，標準誤差は`0.010`。\n",
    "\n",
    "次に，`wage`が`20`以上の観測値を取り除き，`20`未満のサンプルだけで推計する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "htv_20 = htv.query('wage < 20')  # データの抽出\n",
    "\n",
    "print(f'切断前の標本の大きさ：{len(htv)}')\n",
    "print(f'切断前の標本の大きさ：{len(htv_20)}')\n",
    "print(f'削除された標本の大きさ：{len(htv)-len(htv_20)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "164のサンプルが取り除かれた。これにより，ランダムな標本ではなくなっていおり，GM仮定２が満たされていない。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "res_ols_20 = smf.ols(formula_trunc,data=htv_20).fit()\n",
    "\n",
    "print(res_ols_20.summary().tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`educ`の係数は`0.0579`に大きく下落している。切断データをOLS推定すると（ゼロ方向に）バイアスが発生することを示している。\n",
    "\n",
    "---\n",
    "次に，切断回帰推定をおこなう。\n",
    "* 右切断なので`right`に数値を設定する。\n",
    "* 説明変数が対数化されているため，それに合わせて`right=np.log(20)`とする。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "res_trunc = Truncreg.from_formula(formula_trunc,right=np.log(20),data=htv_20).fit()\n",
    "\n",
    "print(res_trunc.summary().tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`educ`の係数は`0.1060`になり，切断される前の標本をOLS推定した際の係数と近い。\n",
    "\n",
    "（コメント）\n",
    "\n",
    "このように切断回帰は，切断データを使い`y`（賃金）に対する`x`（教育）の効果を推定可能とする。一方で，切断されたデータの中での`y`に対する`x`の効果に興味がある場合，バイアスが発生するため，その限界効果の絶対値は$\\left|\\hat{\\beta}_{\\text{Truncreg}}\\right|$よりも低くなる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "＜`Log(Sigma)`について＞\n",
    "* 誤差項は正規分布に従うと仮定され，最尤法により変数の係数$\\beta$と誤差項の標準偏差$\\sigma$が推定される。誤差項の標準偏差の推定値または回帰の標準偏差（`Sigma` = Standard Error of Regression）の対数が`Log(Sigma)`である。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "`dir()`もしくは`py4macro`モジュールの`see()`関数を使い推定結果`res_trunc`の属性とメソッドを確認してみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true,
    "tags": [
     "output_scroll"
    ]
   },
   "outputs": [],
   "source": [
    "py4macro.see(res_trunc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "対数最尤関数の値"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "res_trunc.llf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "疑似決定係数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "res_trunc.prsquared"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "全ての説明変数（定数項以外）が０の場合の推定結果"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "print(res_trunc.result_null.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 検定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "次に検定方法について説明する。\n",
    "\n",
    "**Wald検定の例１**\n",
    "\n",
    "$H_0$：定数項以外の全てのスロープ係数は`0`\n",
    "\n",
    "$H_A$：少なくとも１つの係数は`0`ではない\n",
    "\n",
    "この検定のために結果のメソッド`wald_test()`を使うが説明変数に含まれている`Log(Sigma)`は残る必要がある。従って，まず定数項と`Log(Sigma)`以外の係数名を`slopes_all`に割り当てる。\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "slopes_all = res_trunc.model.exog_names[1:-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "結果`res_trunc`に属性`model`があり，その中に説明変数名の属性`exog_names`を使っている。定数項と最後にくる`Log(Sigma)`を省くために`[1:-1]`を指定している。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Wald検定\n",
    "print( res_trunc.wald_test(slopes_all, scalar=True).summary() )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "$p$値は非常に低いので`1%`の有意水準でも帰無仮説を棄却できる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "２つ目の例を考えよう。\n",
    "\n",
    "**Wald検定の例２**\n",
    "\n",
    "$H_0$：`educ`+`abil`$=$`exper`\n",
    "\n",
    "$H_A$：`educ`+`abil`$\\neq$`exper`\n",
    "\n",
    "次の方法でおこなう。\n",
    "1. 制約式を文字列で設定する：`educ+abil=exper`\n",
    "1. 推定結果のメソッド`wald_test`に制約式を引数として実行する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "print(res_trunc.wald_test('educ+abil=exper', scalar=True).summary()\n",
    "     )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`1%`の有意水準でも帰無仮説を棄却できる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### 予測値と残差"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "次に２つの属性を紹介する。\n",
    "* `.fittedvalues`：以下の式で与えられる線形の予測値\n",
    "\n",
    "    $$\\hat{y}_i=\\hat{\\beta}_0+\\hat{\\beta}_1x_i$$\n",
    "    \n",
    "* `.resid`：以下の式で与えられる線形の残差\n",
    "\n",
    "    $$\\hat{u}_i=y_i-\\hat{y}_i$$\n",
    "    \n",
    "まず予測値の平均・最小値・最大値を計算してみる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "y_hat = res_trunc.fittedvalues\n",
    "\n",
    "print(f'最小値：{y_hat.min()}\\n平均：{y_hat.mean()}\\n最大値：{y_hat.max()}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "次に残差を図示する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "u_hat = res_trunc.resid\n",
    "plt.scatter(y_hat,u_hat)\n",
    "plt.xlabel('y_hat')\n",
    "plt.ylabel('u_hat')\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "データの切断による影響が右上に見て取れる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true,
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "from scipy.stats import truncnorm\n",
    "plt.scatter(htv_20['lwage'],y_hat)\n",
    "plt.ylim(0,10)\n",
    "plt.xlim(0,4)\n",
    "\n",
    "X = res_trunc.fittedvalues\n",
    "s = res_trunc.params[-1]\n",
    "right = (np.log(20) - X)/np.exp(s)\n",
    "imr = truncnorm.logpdf(-X/np.exp(s),a=-np.inf,b=right)\n",
    "yexp = X-np.exp(s)*imr\n",
    "plt.plot(htv_20['lwage'],yexp,'or')\n",
    "# norm.cdf(X, loc=l, scale=s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 誤差項の仮定について"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "切断回帰推定は最尤法を使っており，次の仮定が重要な役割を果たしている。\n",
    "\n",
    "$$u|x\\sim N(0,\\sigma^2)$$\n",
    "\n",
    "即ち，説明変数を所与とすると誤差項は正規分布に従い分散は一定であるという仮定である。\n",
    "\n",
    "正規分布に関しての検定もあるが`py4etrics`には実装されていないので，この仮定の重要性を念頭に置いて推定する必要がある。\n",
    "\n",
    "次に不均一分散について考える。確かめるためには検定をする必要があるが，`py4etrics`には実装されていない。その代わりに「目安」として，通常の係数の標準偏差と不均一分散頑健標準偏差の変化を調べてみる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# 通常の標準偏差\n",
    "tr0 = res_trunc.bse\n",
    "\n",
    "# 不均一分散頑健標準偏差\n",
    "tr1 = Truncreg.from_formula(formula_trunc,\n",
    "                            right=np.log(20),\n",
    "                            data=htv_20).fit(cov_type='HC1',\n",
    "                            disp=False).bse\n",
    "\n",
    "# 不均一分散頑健標準編を使った場合の標準偏差の変化率（％）\n",
    "(100*(tr1-tr0)/tr0)[:-1]   # 最後は Log(Sigma) なので省く"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "標準偏差が減少した変数と増加したし変数がある。特別に大きくないように見えるが，これは目安であることを念頭に置いておく必要がある。\n",
    "\n",
    "（注意）不均一分散の下での最尤推定\n",
    "* 推定量は一致性を満たさない\n",
    "* 標準誤差も一致性を満たさない\n",
    "    * 不均一分散頑健標準誤差を使うことが推奨されることがあるが（研究論文でもそうする研究者も多い），もともと係数の推定量が一致性を満たさないため，`cov_type`で指定する不均一分散頑健標準誤差の有用性に疑問が残る。（[参照](https://davegiles.blogspot.com/2013/05/robust-standard-errors-for-nonlinear.html)）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Tobitモデル"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 説明"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "打ち切りデータを理解するために，切断データと比べて異なる点を説明する。\n",
    "* 切断データではデータセットに$(x_i,y_i),\\;i\\in\\cal{N}$が存在しないが，打ち切りデータには含まれる。しかし，$y_i$が打ち切りの下限や上限の値と等しくなる。\n",
    "\n",
    "例として女性の労働供給を考えよう。働いている女性は多いが労働市場から退出している女性も多いのが実状である。即ち，職を持つ女性の労働供給（例えば，一週間の労働時間）は正の値をとるが，働いていない女性の労働供給は`0`となる。これは数式で次のように表すことができる。\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "&y^{*}=\\beta_0+\\beta_1x+u\\\\\n",
    "    &\\begin{cases}\n",
    "    y=y^{*}&\\quad\\text{if }y^{*}>0\\\\\n",
    "    y=0&\\quad\\text{if }y^{*}\\leq0\n",
    "    \\end{cases}\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "* $y^{*}$：潜在変数（例えば，効用と考えても良い）\n",
    "    * $y^{*}>0$の場合に実際に働き，労働供給は$y=y^{*}$となる。\n",
    "        * $y^{*}$を効用と解釈すると，効用が正の場合，労働者は働く（$y>0$）。その場合，効用は労働時間と同じと仮定する。（$y=y^{*}>0$）。\n",
    "    * $y^{*}\\leq 0$の場合に働かないので$y=0$となる。\n",
    "        * $y^{*}$を効用と解釈すると，効用が`0`以下の場合，労働者は働かない（$y=0$）。その場合，効用と労働時間は異なると仮定する。\n",
    "    * $y^{*}$は観察不可能\n",
    "        * $y^{*}$を効用と解釈すると，最もらいし仮定である。\n",
    "* $x$：労働供給に関する決定要因（例えば，教育年数）\n",
    "* $y$：実際の労働供給（非負の連続変数）\n",
    "* $u$：労働供給に関するランダムな要素（例えば，好み）\n",
    "\n",
    "    $$u|x\\sim \\text{Normal}(0,\\sigma^2)$$\n",
    "    \n",
    "    * この仮定が非常に重要であり，計算の前提となる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "（コメント）\n",
    "* 上の例では，女性が労働市場に参加するかしないかによって，$y$が正の値もしくは`0`を取る。即ち，`0`が下限になっている。上限の例として，人気歌手のコンサート・チケットがあげられる。チケット数は限られており，売り切れた場合の需要は上限を上回る。また，下限・上限の両方がある場合として大学入試が挙げられる。下限はセンター試験などでの「足切り」であり，上限は定員数でる。\n",
    "* 労働供給の例では，女性の選択の結果として$y$の値が観察される。これはミクロ経済学でおなじみの端点解の例である。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "＜結果＞\n",
    "* $y>0$と$y=0$の両方のデータを使ったOLS推定量は不偏性・一致性を満たさない。\n",
    "\n",
    "---\n",
    "以下ではまずこの結果をシミュレーションを使って示し，解決策について説明する。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### シミュレーション"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# データの大きさ \n",
    "n = 100\n",
    "\n",
    "# y*を決定するx\n",
    "x = np.sort(norm.rvs(0,3,size=n))  # ランダム変数を生成し昇順に並べる\n",
    "\n",
    "# y*を生成\n",
    "y_star = x + norm.rvs(0,3,size=n)\n",
    "\n",
    "# yを生成\n",
    "y = y_star.copy()  #  copy()はコピーを作るメソッド\n",
    "y[y_star < 0] = 0  # y_star<0が場合，0を代入する\n",
    "\n",
    "# DataFrame\n",
    "df = pd.DataFrame({'Y':y, 'Y_star':y_star, 'X':x})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "母集団回帰式は次式で与えられる。\n",
    "\n",
    "$$y^{*}=\\beta_0+\\beta_1x+u$$\n",
    "\n",
    "ここで$y^{*}$を効用と解釈すると分かりやすいだろう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "formula_star = 'Y_star ~ X'\n",
    "\n",
    "result_star= smf.ols(formula_star, data=df).fit()\n",
    "\n",
    "b0_star,b1_star = result_star.params\n",
    "\n",
    "print(result_star.params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "$y\\geq 0$を使ったOLSの回帰式は次式となる。\n",
    "\n",
    "$$y=\\beta_0+\\beta_1x+u$$\n",
    "\n",
    "ここで$y$は観察できる労働時間である。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "formula_sample = 'Y ~ X'\n",
    "\n",
    "result_corner = smf.ols(formula_sample, data=df).fit()\n",
    "\n",
    "b0_corner, b1_corner = result_corner.params\n",
    "\n",
    "print(result_corner.params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "母集団回帰式とOLS回帰の推定値は大きく異なる。図示して確かめてみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# y_starの散布図\n",
    "plt.scatter(x, y_star, \n",
    "            facecolors='none', \n",
    "            edgecolors='gray',\n",
    "            label=r'$y^{*}$')\n",
    "\n",
    "# yの散布図\n",
    "plt.scatter(x, y, \n",
    "            facecolors='gray', \n",
    "            label=r'$y$')\n",
    "\n",
    "# 母集団OLS\n",
    "plt.plot(x, b0_star+b1_star*x,\n",
    "         color='k',\n",
    "         linewidth=3,\n",
    "         label='潜在変数（効用）: $y^*$')\n",
    "\n",
    "# y>=0のOLS\n",
    "plt.plot(x, b0_corner+b1_corner*x,\n",
    "         color='r',\n",
    "         lw=3,\n",
    "         label=r'観察可能なサンプル: $y\\geq 0$')\n",
    "\n",
    "plt.xlabel('x')\n",
    "plt.ylabel(r'$y$ and $y^{*}$')\n",
    "plt.legend()\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "この場合，$y\\leq0$の値を全て$y=0$としてOLS推定しているため，不偏性・一致性が満たされないのは直感的に理解できる。実際，上で説明したように標本回帰式は\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "&y^{*}=\\beta_0+\\beta_1x+u\\\\\n",
    "    &\\begin{cases}\n",
    "    y=y^{*}&\\quad\\text{if }y^{*}>0\\\\\n",
    "    y=0&\\quad\\text{if }y^{*}\\leq0\n",
    "    \\end{cases}\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "であるが，$y\\geq 0$の下での$y$の期待値 $\\text{E}(y|y>0,x)$ は複雑な非線形式なり，線形を仮定するOLS推定方で捉えることは不可能である。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### `Tobit`モジュールの使い方"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "残念ながら`statsmodels`と`linearmodels`にも`Tobit`推定のモジュールがない。その代わりに著者が作成した`py4etrics`パッケージの`tobit`モジュールを使い推定する。このモジュールは`statsmodels`の`GenericMaximumLikelihoodModel`を使い実装したものである。使用する上で[このサイト](https://www.statsmodels.org/devel/examples/notebooks/generated/generic_mle.html)にある次の点に注意する必要がある。\n",
    "* `R`の推定値と小数点第４位まで同じになるが，標準偏差は小数点第２位までが同じとなる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "＜使い方＞\n",
    "\n",
    "基本的に`statsmodels`の`ols`と同じだが，追加的な操作とオプションがある。\n",
    "\n",
    "1. 下限・上限の設定：被説明変数`y`の値に従って`Numpy`の`array`もしくは`Pandas`の`Series`を作る。\n",
    "    * 下限がある場合：`-1`\n",
    "    * 上限がある場合：`1`\n",
    "    * 上限と下限の両方がない場合：`0`\n",
    "2. 推定式を決める\n",
    "```\n",
    "    formula = 'y ~ 1 + x'\n",
    "```\n",
    "3. `Tobit`の`from_formula`モジュールを使って推定\n",
    "```\n",
    "    Tobit.from_formula(formula, cens=<A>, left=<B>, right=<C>, data=<D>).fit()\n",
    "```\n",
    "   ここで \n",
    "   \n",
    "    * `cens`：ステップ１で作成した下限・上限の設定が含まれる`array`もしくは`Series`を指定する。\n",
    "    * `left`:下限の値（デフォルトは`0`）\n",
    "        * ステップ１で`-1`が設定されている場合のみ有効となる。\n",
    "    * `right`:上限の値（デフォルトは`0`）\n",
    "        * ステップ１で`1`が設定されている場合のみ有効となる。\n",
    "    * `deta`:データの指定\n",
    "\n",
    "（コメント）\n",
    "*`Logit`や`Probit`と同じように，非線形モデルなため最尤法を使い推定する。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "```{warning}\n",
    "`cens`に指定する`array`の要素数，もしくは`Series`の行数が，`data`で指定する`DataFrame`の行数と同じになる必要がある。異なる場合はエラーが発生する。必ず`len()`関数を使って確認しよう！\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### Tobitモデルの推定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "例として、[](sec:19-case1)で使った`mroz`のデータを使う推定を考えよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "mroz = wooldridge.data('mroz')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "女性の労働供給のモデルを考え，供給量`hours`を被説明変数とする。特に，`hours`は`0`が下限となっているため`Tobit`モデルが妥当だと考えられる。労働時間`hours`を図示してみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "plt.hist(mroz['hours'],bins=20)\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`0`に多くの観測値があることがわかる。\n",
    "\n",
    "まず，下限の値を設定する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "left = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "次に下限を示す`array`を作成する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "cond = (mroz['hours'] == left)  #  フィルターの作成\n",
    "\n",
    "censor = np.zeros((len(mroz)))  # 0のarrayの作成\n",
    "\n",
    "censor[cond] = -1  #  条件に合わせて-1を代入\n",
    "\n",
    "pd.Series(censor).value_counts()  # Serieに変換し，内訳の確認"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "次のコードでも同じ結果を得ることができる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "censor = mroz['hours'].apply(lambda x: -1 if x==left else 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "まず先に`censor`と`mroz`の行数が同じか確認しよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "len(censor) == len(mroz)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "確認できたので，推定してみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "formula = 'hours ~ 1 + nwifeinc + educ + exper + expersq + age + kidslt6 + kidsge6'\n",
    "\n",
    "res_tobit = Tobit.from_formula(formula,\n",
    "                               cens=censor,\n",
    "                               left=0,\n",
    "                               data=mroz).fit()\n",
    "\n",
    "print(res_tobit.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "係数の解釈はOLSと同じようにおこなう。例えば，`educ`の推定値は約`80`なので，教育年数が一年増えると平均で労働時間が年間`80`時間増加することを示している。また`Log(Sigma)`は切断回帰モデルの場合の解釈と同じである。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "`dir()`や`py4macro.see()`を使うことにより属性とメソッドを確認できる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true,
    "tags": [
     "output_scroll"
    ]
   },
   "outputs": [],
   "source": [
    "py4macro.see(res_tobit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "この中にあるメソッドを使い，次の節では検定をおこなう。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 検定と属性"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`Wald検定`を考えよう\n",
    "\n",
    "**Wald検定の例**\n",
    "\n",
    "$H_0$：`exper` $=$ `expersq` $=0$ & `kidslt6`$=$`kidsge6`\n",
    "\n",
    "$H_A$：$H_0$は成立しない\n",
    "\n",
    "検定方法は切断回帰の`Truncreg`モジュールと同じである。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "print( res_tobit.wald_test('exper=expersq=0, kidslt6=kidsge6', scalar=True).summary()\n",
    "     )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "$p$値は非常に低いため，`1%`有意水準でも帰無仮説を棄却できる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "次の３つの属性も有用である。\n",
    "* `.fittedvalues`：以下の式で与えられる線形の予測値 $\\hat{y}^{*}$\n",
    "\n",
    "    $$\\hat{y}_i^{*}=\\hat{\\beta}_0+\\hat{\\beta}_1x_i$$\n",
    "    \n",
    "* `.fitted_endog`：被説明変数の予測値 $\\text{E}(y|x)$ （正規分布に基づいた非線形）\n",
    "* `.resid`：以下の式で与えられる線形の残差\n",
    "\n",
    "    $$\\hat{u}_i=y_i-\\hat{y}_i^{*}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "y_star_hat = res_tobit.fittedvalues\n",
    "\n",
    "u_hat = res_tobit.resid\n",
    "\n",
    "plt.scatter(y_star_hat,u_hat)\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "図の左下は切り打ちデータを反映している。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 残差について"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "\n",
    "切断回帰推定と同じように最尤法を使っているので，推定量の標準偏差の変化を使い残差の不均一について考えてみる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# 通常の標準偏差\n",
    "to0 = res_tobit.bse\n",
    "\n",
    "# 不均一分散標準偏差\n",
    "res_tobit_HC1 = Tobit.from_formula(formula,cens=censor,left=left,\n",
    "                         data=mroz).fit(cov_type='HC1',disp=False)\n",
    "to1 = res_tobit_HC1.bse\n",
    "\n",
    "# 不均一分算標準偏差を使った場合の標準偏差の変化（％）\n",
    "(100*(to1-to0)/to0)[:-1]     # Log(Sigma)を省く"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "全ての係数の標準誤差は大きくは変化しいない。それほど心配する必要もないかも知れない。\n",
    "\n",
    "（注意）不均一分散の下でお最尤推定\n",
    "* 推定量は一致性を満たさない\n",
    "* 標準誤差も一致性を満たさない\n",
    "    * 不均一分散頑健標準誤差を使うことが推奨されることがあるが（研究論文でもそうする研究者も多い），もともと係数の推定量が一致性を満たさないため，`cov_type`で指定する不均一分散頑健標準誤差の有用性に疑問が残る。（[参照](https://davegiles.blogspot.com/2013/05/robust-standard-errors-for-nonlinear.html)）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Heckitモデル"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### 説明"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`Heckit`モデルは切断回帰モデルの拡張版であり，選択バイアス問題が疑われる場合に威力を発揮する。例を使って選択バイアス問題を考えてみよう。\n",
    "\n",
    "＜選択バイアス問題：例１＞\n",
    "\n",
    "大学４年生の学力を測るために，単位取得できる同じ試験を日本全国の大学４年生に受けさせるとしよう。この場合，勉強時間が短い学生は受験しない傾向にあるとしよう。\n",
    "\n",
    "２つのシナリオ\n",
    "* シナリオ１：学生を無作為に選び（強制的に）受験させる。\n",
    "* シナリオ２：希望する学生だけに受けさせる。\n",
    "\n",
    "結果\n",
    "* シナリオ１：ランダム抽出なので平均点は全体像を反映している。\n",
    "* シナリオ２：勉強時間が短い学生は受験しなくなり，平均点が上がることになる。全体像を歪める結果となる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "＜選択バイアス問題：例２＞\n",
    "\n",
    "賃金に対する大学教育１年間の平均的効果を検証するとしよう。サンプルは大学卒業生と大学に進学しなかった高校卒業生。目的は全体像の把握であり，以下を考える。\n",
    "\n",
    "$$\\dfrac{W_{\\text{大}}-W_{\\text{高}}}{4}\\qquad (式１)$$\n",
    "\n",
    "* $W_{\\text{大}}$：大卒の平均賃金\n",
    "* $W_{\\text{高}}$：高卒の平均賃金\n",
    "\n",
    "次の仮定を置く：\n",
    "* 教育が高いほど賃金は高い。\n",
    "* 他の条件が同じであれば，教育が高いほど働く傾向にある（機会費用が高いため）。\n",
    "\n",
    "２つのシナリオ\n",
    "* 起こりえないシナリオ：卒業生を無作為に選び（強制的に）働かせて賃金のデータを集める。\n",
    "* 現実的なシナリオ：自由に労働者が働くか働かないかを決め，働いている労働者のデータを集める。\n",
    "\n",
    "結果：\n",
    "* 起こりえないシナリオ：ランダム抽出なので(式１)は全体像を反映している。\n",
    "* 現実的なシナリオ：教育水準が低い人（高卒）ほど働かな人いが多い傾向にある。特に，賃金が低い人ほど労働市場から退出する場合（労働供給の減少），高卒の平均賃金$W_{\\text{高}}$は上昇することになり，(式１)は下落する。大学教育１年の効果は，低く推定され全体像を歪める結果となる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "これらの例が示唆するように，選択問題を無視して単純にOLS推定しても推定量にバイアスが発生する可能性がある。この問題に対処する推定方法の１つが`Heckit`モデルである。`Heckit`モデルは２段階で推定する。\n",
    "\n",
    "＜第１段階＞\n",
    "\n",
    "`Probit`モデルを使った２項選択に関する推定（例えば，労働市場に参加するしないの選択）\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "&z_i^{*}=\\alpha_0+\\alpha_1w_i+u_i\\qquad\\text{(式５)}\\\\\n",
    "    &\\begin{cases}\n",
    "        z_i=1&\\quad\\text{if }z_i^{*}>0\\\\\n",
    "        z_i=0&\\quad\\text{if }z_i^{*}\\leq0\n",
    "    \\end{cases}\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "* $z_i^{*}$：選択に関する潜在変数（例えば，効用）\n",
    "* $z_i$：選択を示す指標（例えば，１＝働く，０＝働かない）\n",
    "* $w_i$：選択に影響する要因（例えば，働く時間に影響を及ぼす要因として幼児の数）\n",
    "\n",
    "＜第２段階＞\n",
    "\n",
    "第一段階の結果を使い`OLS`推定（例えば，賃金の推定）\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "&y_i^{*}=\\beta_0+\\beta_1x_i+\\rho\\sigma_e\\lambda\\left(\\hat{k}_i\\right)+e_i\\qquad\\text{(式６)}\\\\\n",
    "    &\\begin{cases}\n",
    "        y_i=y_i^{*}&\\quad\\text{if }z_i=1\\\\\n",
    "        y_i\\text{は切断される}&\\quad\\text{if }z_i= 0\n",
    "    \\end{cases}\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "* $y_i^{*}$：興味がある変数（例えば，労働者の賃金）\n",
    "* $y_i$：実際に観測される$y_i^{*}$の値\n",
    "* $x_i$：$y_i^{*}$に影響する要因（例えば，教育，経験）\n",
    "\n",
    "ここで\n",
    "* $\\hat{k}_i$：第一段階における$\\dfrac{A-\\hat{\\alpha}_0-\\hat{\\alpha}_1w_i}{\\sigma_u}$の推定量\n",
    "* $A=0$：下限の値\n",
    "* $\\sigma_e$：$e_i$の標準偏差\n",
    "* $\\sigma_u$：$u_i$の標準偏差\n",
    "* $\\rho=\\text{Cov}(u_i,e_i)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "（コメント）\n",
    "* ある仮定のもとで`Heckit`推定量は一致性を満たす。\n",
    "* `Heckit`を使わずに，(式６)を直接OLS推定すると$\\lambda\\left(\\hat{w}_i\\right)$を無視することになり，欠落変数バイアスが発生する。\n",
    "* $\\rho=0$の場合，(式６)を直接OLS推定しても欠落変数バイアスは発生しない。この場合，`Heckit`モデルを使う必要はない。即ち，$\\rho\\sigma_e$のOLS推定量でバイアスが発生しいるか確認できる。\n",
    "\n",
    "（注意）\n",
    "* 上の説明では，$w_i$も$x_i$も１変数として説明したが，実際には複数の変数を使うことになる。その際，第１段階の説明変数（上の例では，$w_i$）には第２段階の説明変数に**ない**変数を１つ以上入れる必要がある。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### `Heckit`モジュールの使い方"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`statsmodels`も`linearmodels`も`Heckit`推定の正式モジュールがない。その代わり`statsmodels`に正式にmergeされていない`Heckman`モジュールに基づいて著者が作成した`heckit`モジュールを使う。これにより上述の説明したステップに沿って自動で推定可能となる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "＜使い方：ステップ１〜６＞\n",
    "\n",
    "今まで`statsmodels`を使う場合，`from_formula`メソッドを使ったが，それを使わない方法もある。`Heckit`の使い方はその方法に沿っている。\n",
    "\n",
    "1. 変数を準備する前準備\n",
    "    * (式５)と(式６)を使った説明で，第１段階の被説明変数$z_i$と第２段階の被説明変数$y_i$は次のように連動していることに留意しよう。\n",
    "        $$\n",
    "        \\begin{align*}\n",
    "            &z_i=1\\;\\Rightarrow\\; y_i=y^*_i\\\\\n",
    "            &z_i=0\\;\\Rightarrow\\; y_i\\text{は切断}\n",
    "        \\end{align*}\n",
    "        $$\n",
    "    * この点を利用し，`Heckit`モジュールでは第１段階の被説明変数$z_i$を指定する必要がない。その代わりに，$z_i=0$の場合には$y_i=$`NaN`となるように設定する必要がある。\n",
    "1. 第１段階と第２段階で使う全ての変数が入った`DataFrame`を作成する。以下の説明では`df`と呼ぶ。\n",
    "1. 第２段階の被説明変数を`endog`として定義する。例えば，`df`の列`y`にある場合，以下のようにする。\n",
    "    ```\n",
    "        endog = df.loc[:,'y']\n",
    "    ```\n",
    "1. 第２段階の説明変数だけを抜き取って`exog`という`DataFrame`を作成し，それに定数項の列を加える。例えば，`x1`，`x2`，`x3`が該当する変数だとすると，以下のようにする。また`exog`に定数項を加える。\n",
    "    ```\n",
    "        exog = df.loc[:,[x1,x2,x3]]\n",
    "        exog['Intercept'] = 1.0\n",
    "    ```\n",
    "1. 第１段階の説明変数だけを抜き取って`exog_select`という`DataFrame`を作成し，それに定数項の列を加える。例えば，`w1`，`w2`，`w3`が該当する変数だとすると，以下のようにする。また`exog_select`に定数項を加える。\n",
    "    ```\n",
    "        exog_select = df.loc[:,[w1,w2,w3]]\n",
    "        exog_select['Intercept'] = 1.0\n",
    "    ```\n",
    "1. 以下のように`Heckit`を使い推定する。\n",
    "    ```\n",
    "        Heckit(endog, exog, exog_select).fit()\n",
    "    ```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "```{warning}\n",
    "`endog`、`exog`、`exog_select`の行数が全て同じになる必要がある。異なる場合はエラーが発生する。必ず`len()`関数を使って確認しよう！\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 推定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "#### `mroz`を使った推定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "上で使った`mroz`のデータセットを使い推定する。\n",
    "* 第１段階の説明変数：`educ`，`exper`，`expersq`，`nwifeinc`，`age`，`kidslt6`，`kidsge6` \n",
    "* 第２段階の被説明変数：`lwage`\n",
    "* 第２段階の説明変数：`educ`，`exper`，`expersq`\n",
    "\n",
    "（目的）\n",
    "* 教育の収益率の推定。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "ステップ１"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`mroz`にある変数`inlf`は，1975年に既婚女性が労働市場に参加した場合は`1`，参加しなかった場合は`0`になるダミー変数である。この変数を使い，変数作成の前準備として次の２点を確認する。\n",
    "* `inlf`=1の場合，`lwage`は`NaN`ではない。\n",
    "* `inlf`=0の場合，`lwage`は`NaN`である。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "mroz.query('inlf == 1')['lwage'].isna().sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "このコードの`.isna()`は`lwage`が`NaN`であれば`True`を返す（`isnan()`ではないことに注意，また`isnull()`でも可）。その合計`.sum()`が`0`なので，「`inlf`=1の場合，`lwage`は`NaN`ではない」ことが確認できた。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "( ~mroz.query('inlf == 0')['lwage'].isnull() ).sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "このコードの`()`の中を考えよう。上のコードと同じように，`mroz.query('inlf == 0')['lwage'].isnull()`は`lwage`が`NaN`であれば`True`を返すが，その先頭に`~`をつけるとその逆の`False`を返すことになる。`~`は「反転」という意味である。その合計`.sum()`が`0`なので，「`inlf`=0の場合，`lwage`は`NaN`である」ことが確認できた。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "ステップ２\n",
    "\n",
    "`mros`をそのまま使う。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "ステップ３\n",
    "\n",
    "第２段階の被説明変数を作成する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "endog = mroz.loc[:,'lwage']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "ステップ４\n",
    "\n",
    "第２段階の説明変数を作成する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "exog = mroz.loc[:,['educ', 'exper', 'expersq']]\n",
    "exog['Intercept'] = 1.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "ステップ５\n",
    "\n",
    "第１段階の説明変数を作成する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "exog_select = mroz.loc[:,['educ', 'exper', 'expersq','nwifeinc', 'age', 'kidslt6', 'kidsge6', ]]\n",
    "exog_select['Intercept'] = 1.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "次のステップに進む前に、`endog`、`exog`、`exog_select`の行数が同じか確認しよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "len(endog) == len(exog) == len(exog_select)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "ステップ６"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "推定を行う際，`fit()`にオプションを追加し不均一分散頑健標準誤差を指定する。\n",
    "* `cov_type_1`：第１段階推定でのオプション\n",
    "* `cov_type_2`：第２段階推定でのオプション\n",
    "\n",
    "（注意）\n",
    "* オプションを追加しない場合は，`nonrobust`がデフォルトとなる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "res_heckit = Heckit(endog, exog, exog_select).fit(cov_type_2='HC1')\n",
    "\n",
    "print(res_heckit.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "* 上段の表：第２段階推定\n",
    "* 中断の表：第１段階推定\n",
    "* 下段の表\n",
    "    * `IMR`：逆ミルズ比（(式６)の$\\lambda\\left(\\hat{w}_i\\right)$）\n",
    "    * `rho`：第１・第２段階の誤差項の共分散（(式６)の$\\rho$）\n",
    "    * `sigma`：第２段階誤差項の標準偏差（(式６)の$\\sigma_e$）\n",
    "\n",
    "（注意）表には通常の標準誤差が表示されている。不均一分散頑健標準誤差は以下で手動で計算する。\n",
    "\n",
    "---\n",
    "第２段階結果の属性とメソッドは`dir()`や`py4macro.see()`で確認できる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "py4macro.see(res_heckit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "例えば，`predict()`は予測値を返す。ただ，この中にはまだ実装されていないものも含まれている。\n",
    "\n",
    "また，この中に`select_res`とあるが，この属性を使い第１段階推定のに関する属性・メソッドを確認できる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "py4macro.see(res_heckit.select_res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "例えば，`fittedvalues`は予測値を返す。次のコードでは基本統計量を表示できる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "print(res_heckit.select_res.summary().tables[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "この表にある`Dep. Variable: y`の`y`は第１段階の被説明変数を表しているが，第２段階の被説明変数`lwage`を使っているためこのような表記になっている。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### 結果の解釈"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "* 第２段階推定の`educ`の係数は`0.1091`であり統計的に有意。\n",
    "\n",
    "この結果を単純なOLSと比べよう。（`lwage`にある`NaN`は自動的に除外される。）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "formula = 'lwage ~ educ + exper + expersq'\n",
    "\n",
    "res = smf.ols(formula, data=mroz).fit(cov_type='HC1')\n",
    "\n",
    "print(res.summary().tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`educ`のOLS推定値は`0.1075`で`Heckit`推定値と大きく変わらない。これは選択バイアスが殆どないことを示唆している。実際，`IMR`（逆ミルズ比）のp値は`0.809`であり，係数は０とする通常の有意水準で帰無仮説を棄却できない。即ち，単純なOLS推定では逆ミルズ比の欠落変数バイアスが発生していないことを意味する。\n",
    "\n",
    "次に切断回帰推定を考えてみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "plt.hist(mroz['wage'].dropna(),bins=20)\n",
    "pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "thresh = np.log(mroz['wage'].min()*0.5)  # 左切断の下限\n",
    "\n",
    "formula = 'lwage ~ 1 + educ + exper + expersq'\n",
    "\n",
    "res_trunc = Truncreg.from_formula(formula, left=thresh,\n",
    "                                  data=mroz.dropna(subset=['lwage'])).fit()\n",
    "\n",
    "print(res.summary().tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "この推定では，`wage`の最小値の50%の値の対数を下限の値に設定している。`wage=0`を下限にしてもよいが，その場合，`np.log(0)`はマイナス無限になり，通常の最尤推定と同じになる。切断回帰推定を使っても結果は変わらない。即ち，選択バイアスが大きな問題ではないことを示唆している。\n",
    "\n",
    "---\n",
    "第１段階推定では`Probit`モデルを使っているが，以下では不均一分散に関して検定を行う。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "het_test_probit(res_heckit.select_res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "帰無仮説は棄却できない。推定量の一致性は満たされていると考えて良いだろう。"
   ]
  }
 ],
 "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"
  },
  "nteract": {
   "version": "0.23.3"
  },
  "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
}