8_Simple_Regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(chap:8-simple)=\n",
    "# 単回帰分析"
   ]
  },
  {
   "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 py4macro    # 属性を調べるsee()関数が含まれる\n",
    "import statsmodels.formula.api as smf  # 回帰分析のパッケージ\n",
    "import wooldridge  # wooldridgeのデータ・パッケージ\n",
    "\n",
    "from scipy.stats import norm    # 正規分布関数\n",
    "from numba import njit  # シミュレーションの計算を高速化する\n",
    "\n",
    "# 警告メッセージを非表示\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## 説明"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "説明変数が１つの単回帰（Simple Regression）を考える。\n",
    "\n",
    "$$y_i=\\beta_0+\\beta_1 x_i + u_i$$\n",
    "\n",
    "ここで使う変数の呼び名は複数あるので混同しないように。\n",
    "* $y_i$： 被説明変数（explained variable）\n",
    "    * 従属変数（dependent variable）\n",
    "    * regressand\n",
    "* $x_i$： 説明変数（explanatory variable）\n",
    "    * 独立変数（independent variable）\n",
    "    * regressor\n",
    "* $u_i$： 誤差項（errors)\n",
    "    * 観測不可能なランダムな要素\n",
    "    * 平均はゼロと仮定する（重回帰分析の章でより詳しく説明する）\n",
    "* $i=1,2,...n$：観測値のインデックス\n",
    "* $\\beta_0$：定数項（観測不可能なため推定する）\n",
    "* $\\beta_1$：スロープ係数（観測不可能なため推定する）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "**最小二乗法推定法（Ordinary Least Squares; OLS）**\n",
    "\n",
    "$$\\sum_{i=1}^n\\left(y_i-\\beta_0-\\beta_1 x_i\\right)^2$$\n",
    "\n",
    "を最小化する$\\beta_0$と$\\beta_1$を選択すると以下を得る：\n",
    "\n",
    "$$\n",
    "\\hat{\\beta}_0=\\bar{y}-\\hat{\\beta}_1\\bar{x},\\qquad\\qquad\n",
    "\\hat{\\beta}_1=\\frac{\\text{Cov}(x,y)}{\\text{Var}(x)}\n",
    "$$\n",
    "\n",
    "* `^`（ハット）がある変数（$\\hat{\\beta}_0$と$\\hat{\\beta}_1$）を**OLS推定量**と呼ぶ\n",
    "* `-`（バー）がある変数（例えば，$\\bar{y}$）は平均値\n",
    "* $\\text{Cov}(x,y)=$ $x$と$y$の共分散\n",
    "* $\\text{Var}(x)=$ $x$の分散\n",
    "* $x$のように添え字$i$を省いた場合は，$x_1$から$x_n$までの全ての観測値という意味（$y$も同じ）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "また以下を定義する。\n",
    "* 予測値 $\\hat{y}_i$\n",
    "\n",
    "$$\\hat{y}_i=\\hat{\\beta}_0 + \\hat{\\beta}_1x_i$$\n",
    "\n",
    "* 残差 $\\hat{u}_i$\n",
    "\n",
    "$$\\hat{u}_i=y_i-\\hat{y}_i$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## データ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Wooldridge (2019)で使われるデータセットは全て`wooldridge`のパッケージに含まれている。そのパッケージからからデータを読み込んで使うことにする。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "````{note}\n",
    "MacではTerminal、WindowsではGit Bashを使い、次のコマンドで`wooldridge`パッケージをインストールできる。\n",
    "```\n",
    "pip install wooldridge\n",
    "```\n",
    "````"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "使う関数は`data()`であり，引数なしで実行するとデータセットのリストが表示される。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "wooldridge.data()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "この中にある`wage1`を使う。データ・セットを読み込むには、引数にデータ・セット名を文字列として入力する。返り値は`DataFrame`であり、\n",
    "回帰分析に使うために`df`に割り当てる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "df = wooldridge.data('wage1')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "オプション`description=True`を追加するとデータの説明を表示することができる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "wooldridge.data('wage1',description=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "このデータを使い次の問題を考える。\n",
    "\n",
    "> 教育水準は賃金にどのような影響を与えるのか\n",
    "\n",
    "* 被説明変数：`wage`（平均時給、単位：米ドル）\n",
    "* 説明変数：`educ`（教育年数、単位：年）\n",
    "\n",
    "まず関係する変数だけを取り出してみる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "df = df.loc[:,['wage','educ']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "最初の3行の表示をしてみよう。\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "df.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "次に`DataFrame`の特徴を確認する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "df.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "また記述統計も必ず確認するように。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "df.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`.head()`，`.tail()`，`.info()`そして`.describe()`を使い`df`の情報を確認し，データのイメージ作りをすることは回帰分析に役立つだろう。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## 「手計算」"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`Python`では回帰分析の計算を自動化するパッケージが提供されており、後ほどその１つを紹介する。その前に、自動計算にはどのようなコードが使われているかを感じてもらうために「手計算」で、即ち、上で示した$\\hat{\\beta}_0$と$\\hat{\\beta}_1$の式を使い係数の推定値を計算してみる。まず分散共分散行列を計算する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "mat_wage_educ = np.cov(df['wage'],df['educ'])\n",
    "mat_wage_educ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "共分散の値を`cov_wage_edu`に割り当てる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "cov_wage_educ = mat_wage_educ[0, 1]\n",
    "cov_wage_educ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`wage`と`educ`の分散は次のコードで取得できる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "mat_wage_educ[0,0], mat_wage_educ[1,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "次のように計算しても同じ値となる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "var_wage = df['wage'].var(ddof=1)\n",
    "var_educ = df['educ'].var(ddof=1)\n",
    "var_wage, var_educ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "それぞれの変数の平均を計算する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "mean_wage = df['wage'].mean()\n",
    "mean_educ = df['educ'].mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "上の式に従って$\\beta_1$の推定値 $\\left(\\hat{\\beta}_1\\right)$を計算する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "b1hat = cov_wage_educ / var_educ\n",
    "b1hat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "**＜結果の解釈＞**\n",
    "\n",
    "`b1hat`は次式の推定結果である。\n",
    "\n",
    "$$y = \\beta_0+\\beta_1x+u$$\n",
    "\n",
    "この式を$x$について微分すると\n",
    "\n",
    "$$\n",
    "\\dfrac{dy}{dx}=\\dfrac{y\\text{の変化}}{x\\text{の増加}}\n",
    "=\n",
    "\\left(\n",
    "\\begin{array}{l}\n",
    "    x\\text{が1単位増加した}\\\\\n",
    "    \\text{場合の}yの\\text{変化}\n",
    "    \\phantom{春}\n",
    "\\end{array}\n",
    "\\right)\n",
    "=\\beta_1\n",
    "$$\n",
    "\n",
    "となる。即ち，$\\beta_1$は教育年数が１年増加した場合，平均時給は何ドル変化するかを示している。OLS推定結果は教育年数が１年増えると賃金は約0.54ドル上昇することを示している。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "同様に$\\beta_0$の推定値 $\\left(\\hat{\\beta}_0\\right)$も計算する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "b0hat = mean_wage - b1hat * mean_educ\n",
    "b0hat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "source": [
    "## `statsmodels`を使う"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(sec:8-statsmodels)=\n",
    "## `statsmodels`を使う"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Python`では簡単なコードで回帰分析に関する多くの計算を行うことができるパッケージが提供されている。その内、代表的なパッケージが`statsmodels`である。次の順番で使い方を説明する。\n",
    "* `ols`サブパッケージを読み込む\n",
    "* 回帰式を定義する\n",
    "* 推定値の自動計算"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 1: サブパッケージの読み込み**\n",
    "\n",
    "`statsmodels`には多くのサブパッケージや関数があるが，ここでは回帰式を文字列で書くことができる`ols`関数をつかう。（他のコードの書き方もあるが，他のトピックで扱う。）`ols`はこのページの冒頭で次のコードでインポートされている。\n",
    "```\n",
    "from statsmodels.formula.api import ols\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 2: 回帰式の定義**\n",
    "\n",
    "* 式を`''`で囲む（文字列）\n",
    "* 定数項である（$\\beta_0$）は自動的に挿入される。\n",
    "    * 定数を省きたい場合は式に`-1`を加える。\n",
    "    * 後に扱う`linearmodels`パッケージでは，定数が自動的に挿入されないので要注意。\n",
    "\n",
    "---\n",
    "定数項がある場合：\n",
    "```\n",
    "被説明変数 ~ 説明変数\n",
    "```\n",
    "定数項がない場合：\n",
    "```\n",
    "被説明変数 ~ -1 + 説明変数\n",
    "```\n",
    "\n",
    "---\n",
    "`wage`を被説明変数にし，定数項を加えて`educ`を説明変数をする場合，次のように書くことになる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "formula = 'wage ~ educ'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 3: 自動計算**\n",
    "\n",
    "まず自動計算の準備として計算の対象となるもの（インスタンスと呼ばれるオブジェクト）を生成し、ここでは\"model\"の省略形の変数`mod`に割り当てる（`mod`でなくとも好きな変数名でOK）。ここで第１引数は上で定義した回帰式，第２引数`data`は使用するデータをを指定する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod = smf.ols(formula, data=df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`dir()`もしくは`py4macro`モジュールに含まれる`see()`関数を使って，`mod`の属性を調べてみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "py4macro.see(mod)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "この中にOLS関連の計算を一瞬で実行するメソッド`fit()`が準備されており、それを使い結果を`res`に割り当てる（`res`はresultの略）。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = mod.fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上の２つを分けずに一行で書いても良い。\n",
    "```\n",
    "res = ols(formula, data=df).fit()\n",
    "```\n",
    "`res`の属性を確認してみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "py4macro.see(res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "例えば，上のリストの中に`.params`とあるが，この属性を使い係数の推定値を取り出すことができる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "「手計算」での計算結果と同じだと確認できる。\n",
    "\n",
    "また`wage`の予測値は`res.fittedvalues`もしくは`res.predict()`で得ることができる。`res.fittedvalues`は属性であり`res.predict()`はメソッドであるが、後者は引数がない場合は前者と同じ結果を返すことになる。\n",
    "\n",
    "上のリストが示すように`res`には多くの情報が詰まっており，これらを駆使して図を描いたり推定値の検定をおこなうことになる。\n",
    "\n",
    "この中にあるメソッド`summary()`を使うと，主な結果がまとめられた表を出力することができる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(res.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ここで`print()`関数を使ったが、使わなくとも同じ内容を表示できるので試してみよう。\n",
    "\n",
    "また，この表は３つの表から構成されており，必要な部分だけを表示することも可能である。\n",
    "```\n",
    "res.summary().tables[0]\n",
    "res.summary().tables[1]\n",
    "res.summary().tables[2]\n",
    "```\n",
    "例えば，係数の推定値に関する基本的な統計量だけを表示するには以下のようにすれば良い。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(res.summary().tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "````{note}\n",
    "`statsmodels 13.0`では`.summary()`に引数`slim`（デフォルトは`False`）が追加されており，`.tables`を使わずに次のコードで簡略化された表を表示できる。\n",
    "```\n",
    "result.summary(slim=True)\n",
    "```\n",
    "````"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**＜表の説明＞**\n",
    "\n",
    "表は3つのセクションから構成されている。\n",
    "* 上段にはOLS推定の基本的な情報が表示されている。\n",
    "    * 左側\n",
    "        * `Dep. Variable`：被説明変数\n",
    "        * `Model`：モデル\n",
    "        * `Method`：手法\n",
    "        * `Data`：日にち\n",
    "        * `Time`：時間\n",
    "        * `No. Observation`：標本の大きさ\n",
    "        * `Df Residuals`：残差の自由度\n",
    "        * `Df Model`：モデルの自由度（定数項以外の説明変数の数）\n",
    "        * `Covariance Type`：共分散のタイプ\n",
    "    * 右側\n",
    "        * `R-squared`：決定係数\n",
    "        * `adj. R-squared`：自由度調整済み決定係数\n",
    "        * `F-statistic`：$F$統計量\n",
    "        * `Prob (F-statistic)`：$F$値\n",
    "        * `Log-Likelihood`：対数尤度\n",
    "        * `AIC`：赤池情報量規準\n",
    "        * `BIC`：ベイズ情報量規準\n",
    "* 中段には主な推定結果が表示される。\n",
    "    * 列ラベル\n",
    "        * `coef`：係数\n",
    "        * `std err`：標準誤差\n",
    "        * `t`：$t$値\n",
    "        * `P>|t|`：$p$値\n",
    "        * `[0.025,0.975]`：信頼区間（5%)\n",
    "    * 行ラベル\n",
    "        * `Intercept`：定数項\n",
    "        * `tfp_relative`：説明変数（選択する変数によって変わる）\n",
    "* 下段には様々な検定などに関する数値が並んでいる。\n",
    "    * 左側\n",
    "        * `Omnibus`：オムニバス検定統計量（帰無仮説：残差は正規分布に従う）\n",
    "        * `Prob(Omnibus)`：オムニバス検定$p$値（帰無仮説：残差は正規分布に従う）\n",
    "        * `Skew`：残差の歪度（正規分布であれば`0`）\n",
    "        * `Kurtosis`：残差の尖度（正規分布であれば`3`）\n",
    "    * 右側\n",
    "        * `Durbin-Watson`：ダービン・ワトソン統計量（残差の自己相関の検定）\n",
    "        * `Jarque-Bera (JB)`：ジャーク・ベラ検定統計量（帰無仮説：残差は正規分布に従う）\n",
    "        * `Prob(JB)`：ジャーク・ベラ検定$p$値（帰無仮説：残差は正規分布に従う）\n",
    "        * `Cond. No.`：条件指数（Condition Index）の最大値（多重共線性を確認するための指標であり，単回帰分析では無視して良い。一方，`statsmodels`では該当する行列を標準化して計算していないため変数の値の大きさに依存することになり，使い難い指標となっている。重回帰分析においての多重共線性の確認については[ここで](sec:9-vif)説明する手法を使うことを勧める。）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OLS推定結果を割り当てた`res`には，ここで説明した以外に数多くの属性が存在する。それらに関する説明は[このサイト](https://www.statsmodels.org/stable/generated/statsmodels.regression.linear_model.RegressionResults.html)を参照すること。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## 図"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "回帰分析の結果を図示するために、まず予測値の取得する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "wage_hat = res.fittedvalues"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "データの散布図と回帰直線を図示する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "wage_hat = res.fittedvalues\n",
    "plt.scatter('educ', 'wage', data=df)\n",
    "plt.plot(df['educ'],wage_hat, 'r')\n",
    "plt.xlabel('Education')\n",
    "plt.ylabel('Wage')\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## モデルの適合度"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "モデルの適合度を確認するために次を定義する。\n",
    "\n",
    "* SST（Total Sum of Squares; 全変動の平方和）\n",
    "\n",
    "$$\\text{SST}=\\sum_{i=1}^n\\left(y_i-\\bar{y}\\right)^2=(n-1)\\cdot\\text{Var}(y)$$\n",
    "\n",
    "* SSE（Explained Sum of Squares; 回帰変動の平方和）\n",
    "\n",
    "$$\n",
    "\\text{SSE}=\\sum_{i=1}^n\\left(\\hat{y}_i-\\bar{y}\\right)^2=(n-1)\\cdot\\text{Var}\\left(\\hat{y}\\right)\n",
    "$$\n",
    "\n",
    "* SSR（Residual Sum of Squares; 残差変動の平方和）\n",
    "\n",
    "$$\n",
    "\\text{SSR}=\\sum_{i=1}^n\\left(\\hat{u}_i-0\\right)^2=(n-1)\\cdot\\text{Var}\\left(\\hat{u}\\right)\n",
    "$$\n",
    "\n",
    "ここで，$\\text{Var}(z)=\\dfrac{1}{n-1}\\sum_{i=1}^n\\left(z_i-\\bar{z}\\right)^2$。これらを使い決定係数と自由度修正後の決定係数を定義する。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "* 決定係数\n",
    "\n",
    "$$R^2=\\frac{SSE}{SST}=1-\\frac{SSR}{SST}$$\n",
    "\n",
    "$$\n",
    "R^2\n",
    "=\\frac{\\text{Var}\\left(\\hat{y}\\right)}{\\text{Var}\\left(y\\right)}\n",
    "=1-\\frac{\\text{Var}\\left(\\hat{u}\\right)}{\\text{Var}\\left(y\\right)}\n",
    "$$\n",
    "\n",
    "* 自由度修正後の決定係数\n",
    "\n",
    "$$\\overline{R}^2=1-(1-R^2)\\frac{n-1}{n-k-1}$$\n",
    "\n",
    "ここで，$k=$ **定数項以外の説明変数の数**　（上の例では$k=1$）\n",
    "\n",
    "ここで説明した変数の値は，上の例の`res`の属性として取得できる。\n",
    "```\n",
    "n = res.nobs\n",
    "k = res.df_model\n",
    "n-k-1 = res.df_resid\n",
    "SST = res.centered_tss\n",
    "SSE = res.ess\n",
    "SSR = res.ssr\n",
    "決定係数 = res.rsquared\n",
    "修正済み決定係数 = res.rsquared_adj\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "n = res.nobs\n",
    "k = res.df_model\n",
    "# n-k-1 = res.df_resid\n",
    "SST = res.centered_tss\n",
    "SSE = res.ess\n",
    "SSR = res.ssr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "**決定係数**\n",
    "\n",
    "上の定義を使うと次のように計算できる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "1-SSR/SST"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "また次のように計算することも可能である。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "u_hat = res.resid  # 残差\n",
    "wage = df['wage']  # 被説明変数\n",
    "R2 = 1 - np.var(u_hat) / np.var(wage)\n",
    "R2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "しかし「手計算」よりも`res`の属性`.rsquared`を使う方がより簡単である。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "res.rsquared"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "**修正済み決定係数**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "上の定義を使うと次のように計算できる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "adjR2 = 1-(1-R2)*(n-1)/(n-k-1)\n",
    "adjR2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "一方で修正済み決定係数を取得するメソッドは`.rsquared_adj`である。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "res.rsquared_adj"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## 賃金の対数化"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "上の回帰分析では被説明変数`wage`をそのまま使ったが，労働経済学では賃金の対数をとり推定するのが一般的である。そのような式を賃金方程式と呼ぶ。その考えに沿って，ここでは`wage`の自然対数変換をおこない回帰分析をおこなう。まず対数化賃金お変数を作成する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "df['wage_log'] = np.log(df['wage'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "formula = 'wage_log ~ educ'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "res = smf.ols(formula, data=df).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "print(res.summary().tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "**＜結果の解釈＞**\n",
    "\n",
    "`res`は次式の推定結果である。\n",
    "\n",
    "$$\\ln y = \\beta_0+\\beta_1x+u$$\n",
    "\n",
    "この式を$x$について微分すると\n",
    "\n",
    "$$\\dfrac{d\\ln y}{dx}=\\dfrac{dy/y}{dx}=\\dfrac{y\\text{の%変化}}{dx}=\\beta_1$$\n",
    "\n",
    "となる。即ち，$\\beta_1$は$x$（教育年数）が一単位増加した場合，$y$（賃金）が何％変化するかを示している。この値は教育の収益率である。OLS推定結果は教育年数が１年増えると賃金は約8.3%増加することを示している。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "```{note}\n",
    "上の例では`wage_log`の変数を作成し計算したが、変数変換のコードを回帰式に直接書く方法もある。次章で説明することにする。\n",
    "```\n",
    "\n",
    "結果を図示しよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# 予測値の取得\n",
    "wage_log_hat = res.fittedvalues\n",
    "\n",
    "# 散布図と回帰曲線\n",
    "plt.scatter('educ', 'wage_log', data=df)\n",
    "plt.plot(df['educ'], wage_log_hat, 'r')\n",
    "plt.xlabel('Education')\n",
    "plt.ylabel('Log of Wage')\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## シミュレーション：$OLS$推定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "上ではデータを使い「手計算」と`statsmodels`を使い係数を計算したが，ここではOLS推定の「意味」をシミュレーションを使って考えてみることにする。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 母集団の考え方"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "母集団とは分析の対象となる共通の属性などを持つ個人・企業・産業・地域・国などの集合体といえる。ここで重要な点は母集団の大きさとランダム抽出の関係を理解することである。次の２つの例を使って考えてみよう。\n",
    "\n",
    "* 「今年の人間の平均身長」を知りたいとしよう。この場合の母集団は地球上の全ての人間である。しかし，全員の身長を測るのはあまりにも費用が高いため，標本を集め，その平均から「人間の平均身長」を推定することが可能である（[参照](https://ja.wikipedia.org/wiki/大数の法則)）。この場合「標本の大きさ$<$母集団の大きさ」となる。\n",
    "* 「今年のOECD諸国のGDP」についての何らかの統計量を計算したいとしよう。2021年現在でOECDには38カ国のメンバー国があり，GDPのデータは簡単に入手可能である。これは「標本の大きさ$=$母集団の大きさ」という意味なのだろうか。そうだとすれば，標本ではなく母集団全てのデータを扱うことになり，標本の概念である標準誤差などは意味がなくなることになる。しかし「標本の大きさ$<$母集団」と解釈できる２つの考え方がある。\n",
    "    1. 時間は「無限」と考えることができる。この考えに従うと、38カ国のGDPデータも「無限」に存在することになる（OECDが存在する限りという意味での「無限」）。従って、ある時点の38カ国のGDPのデータは「無限」の母集団から抽出されたランダム変数と考えることができる。国と時間の概念で考えると母集団の大きさは非常に大きいと解釈できる。\n",
    "    2. 確率変数とは実現可能な数値から１つを実現させる関数である。またGDP自体は確率変数として解釈する方が自然である（[「重回帰分析」](sec:9-assumptions))で説明する仮定２を参照）。従って，38カ国の**実現可能なGDPの値が母集団**となる。今年度のみでの実現可能なGDPは無数にあるため，この解釈によると「標本の大きさ$<$母集団の大きさ」となる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### プラン"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "**＜Step 1＞**\n",
    "* 母集団の大きさは指定しない。\n",
    "* 母集団のデータを生成する回帰式（母集団回帰式）を定義する。\n",
    "\n",
    "    $$Y = \\beta_0 + \\beta_1 X + U$$\n",
    "\n",
    "    * 説明変数$X$と誤差項$U$はランダム変数であり，従って$Y$もランダム変数である。\n",
    "* 母集団から標本を1組抽出し，標本の大きさを$n$とする。\n",
    "    * 抽出された値はランダム変数の実現値である。\n",
    "    \n",
    "        $$(y_1,x_1),\\;(y_2,x_2),.....(y_n,x_n)$$\n",
    "        \n",
    "    * 仮定\n",
    "        * $X$は正規分布に従う $\\Rightarrow$ $x_1,\\;x_2,..x_n$は$X$の実現値。\n",
    "        * $U$は正規分布に従う $\\Rightarrow$ $u_1,\\;u_2,..u_n$は$U$の実現値。\n",
    "        * $X$と$U$は無相関。\n",
    "        * 実現値$y_1,\\;y_2,..y_n$は母集団回帰式に従って生成される。\n",
    "        \n",
    "            $$y_i = \\beta_0 + \\beta_1 x_i + u_i,\\qquad i=1,2,3,....n$$\n",
    "            \n",
    "    * $\\beta_0$と$\\beta_1$は未知であり，$u_i$も観測不可能\n",
    "\n",
    "**＜Step 2＞**\n",
    "* Step 1で抽出した標本を使い回帰分析を行い，母集団の回帰式のパラメータを推定する。\n",
    "    * $\\hat{\\beta}_0$，$\\hat{\\beta}_1$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### シミュレーション"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "母集団のパラメータを設定する。これはパラメータの真の値であり、この値を推定するのが回帰分析となる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "b0 = 1.0  # 定数項\n",
    "b1 = 1.0  # 説明変数の係数\n",
    "su = 1.0  # 誤差項の標準偏差"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "標本の大きさ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "n = 30"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "標本を抽出するが、次の無作為なプロセスを通して生成されるとする。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "x = norm.rvs(4, 1, size=n)\n",
    "u = norm.rvs(0, su, size=n)\n",
    "y = b0 + b1 * x + u"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "このデータから`DataFrame`を作成する。`()`の中の`{}`は辞書である。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "df_sim = pd.DataFrame({'X':x, 'Y':y})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`statsmodels`を使って`OLS`推定する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "formula = 'Y ~ X'\n",
    "res_sim = smf.ols(formula, data=df_sim).fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "係数の推定値を表示する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "res_sim.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "母集団パラメータと比べると，推定誤差が発生していることがわかる。\n",
    "\n",
    "```{note}\n",
    "ランダム変数を発生させ回帰分析を行っていることに注意しよう。ここで生成したランダム変数は、読者が生成したランダム変数とは異なる。従って，係数の推定値も異なることになる。\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 関数化"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "次に，上で行ったシミュレーションのコードを関数にまとめることにする。\n",
    "\n",
    "`def`を使い`my_sim`という名前の関数を作成する。\n",
    "* 引数\n",
    "    * 標本の大きさ：`n`\n",
    "    * 真の定数項：`b0`（デフォルトは`1`）\n",
    "    * 真のスロープ係数：`b1`（デフォルトは`1`）\n",
    "    * 誤差項の標準偏差：`su`（デフォルトは`1`）\n",
    "* 戻り値\n",
    "    * 定数項とスロープ係数の推定値（`pandas.Series`）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def my_sim(n,b0=1,b1=1,su=1):  # nを指定する，　他はデフォルトを設定\n",
    "    \n",
    "    x = norm.rvs(4, 1, size=n)\n",
    "    u = norm.rvs(0, su, size=n)\n",
    "    y = b0 + b1 * x + u\n",
    "    \n",
    "    df_sim = pd.DataFrame({'X':x, 'Y':y})\n",
    "    \n",
    "    formula = 'Y ~ X'\n",
    "    res_sim = smf.ols(formula, data=df_sim).fit()\n",
    "    \n",
    "    return res_sim.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`my_sim()`関数の中身は，上のコードと殆ど変わらないことに気付くはずである。\n",
    "\n",
    "標本の大きさを指定して関数を実行してみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "my_sim(100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`b0`，`b1`，`su`変えて計算してみよう。\n",
    "\n",
    "---\n",
    "シミュレーションで同じ計算を数十万回おこなう場合，１回のシミュレーションに時間がかかり，それをどのように高速化するかが問題になる。その際に役に立つパッケージが`Numba`であり，１行のコードを追加するだけで計算を驚く程高速化できる。しかし`Numba`には「相性」の問題がある。例えば，上の関数`my_sim()`で使った`statsmodels`，`pandas`，`scipy.stats`は`Numba`と相性が悪いが，`NumPy`は非常に相性が良い。また，`Numba`を使う際，リストよりも`array`を使う方が少しでも高速化を図ることができることも覚えておくと良いだろう。以下では，`Numba`を使う準備として，`NumPy`を使う「手計算」関数を考えることにする。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def my_sim_manual(n,b0=1,b1=1,su=1):  # nは標本の大きさ，　他はデフォルトを設定\n",
    "    \n",
    "    x = np.random.normal(loc=4, scale=1, size=n)   # 以下を参照\n",
    "    u = np.random.normal(loc=0, scale=su, size=n)  # 以下を参照\n",
    "    y = b0 + b1 * x + u\n",
    "    \n",
    "    var_cov = np.cov(y,x)   # 分散共分散行列\n",
    "    cov_yx = var_cov[0,1]   # 共分散\n",
    "    var_y = var_cov[0,0]    # yの分散\n",
    "    var_x = var_cov[1,1]    # xの分散        \n",
    "    avr_y = y.mean()        # yの平均\n",
    "    avr_x = x.mean()        # xの平均\n",
    "    \n",
    "    b1hat = cov_yx / var_x  # スロープ係数の推定値\n",
    "    b0hat = avr_y - b1hat*avr_x  # 定数項の推定値\n",
    "    \n",
    "    return b0hat, b1hat     # 推定値"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "ここで使った`np.random.normal()`について簡単に説明する。`NumPy`にはランダム変数を生成する`random`というサブパッケージがあり，その中の`normal()`を使う。これは`scipy.stats`の`norm.rvs()`と同じ引数を取り，同じ役割を果たす。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "my_sim_manual(100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`my_sim()`と`my_sim_manual()`の推定値が異なるのはランダム変数（`y`，`x`）の値が異なるためであり，それが推定誤差による影響である。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### ループを使ってシミュレーション"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "標本による推定は誤差が発生する。その誤差はどれだけ大きく（小さく）何に依存するかを確認するために`for`ループを使い次のシミュレーションを行うことにする。\n",
    "\n",
    "* 標本の大きさ`n`のデータを`N`組用意し（`N`は標本数），`OLS`推定を`N`回おこなう。\n",
    "* シミュレーション毎に計算した回帰係数の推定値の分布を図示する。\n",
    "\n",
    "まず`def`を使いシミュレーションの関数を作成するが，パラメータの推定値の計算には`statsmodels`を使わずに，`numpy`の関数を使い「手計算」とする。これにより計算を高速化するパッケージ`Numba`を使うことが可能となる。使い方は簡単でデコレーターと呼ばれる`@njit`（又は`@jit`）を関数の上に加えるだけである。これだけで計算速度が数十倍早くなる場合もある。そのために、まず上で使った`my_sim_manual`に`@njit`を加えて実行する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "@njit\n",
    "def my_sim_manual(n,b0,b1,su):\n",
    "    \n",
    "    x = np.random.normal(loc=4, scale=1, size=n)\n",
    "    u = np.random.normal(loc=0, scale=su, size=n)\n",
    "    y = b0 + b1 * x + u\n",
    "    \n",
    "    var_cov = np.cov(y,x)\n",
    "    cov_yx = var_cov[0,1]\n",
    "    var_y = var_cov[0,0]\n",
    "    var_x = var_cov[1,1]\n",
    "    avr_y = y.mean()\n",
    "    avr_x = x.mean()\n",
    "    \n",
    "    b1hat = cov_yx / var_x\n",
    "    b0hat = avr_y - b1hat*avr_x\n",
    "    \n",
    "    return b0hat, b1hat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "この`my_sim_manual()`の内容は，上で作成した`my_sim_manual()`と全く同じである。単に`@njit`の一行だけが付け加えられているだけだが，これにより全く別の高速化された関数に「変身」することになる。\n",
    "\n",
    "次に`my_sim_manual()`を`for`ループで計算するシミュレーション用の関数を用意する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "@njit\n",
    "def my_sim_loop(n,N,b0=1.0,b1=1.0,su=1.0):  # nとN以外はデフォルトを指定\n",
    "    \n",
    "    b0hat_arr = np.zeros(N)  # OLS推定量を入れるarray（0がN個並んでいる）\n",
    "    b1hat_arr = np.zeros(N)  # OLS推定量を入れるarray（0がN個並んでいる）\n",
    "\n",
    "    for i in range(N):  # N 回のループ\n",
    "        \n",
    "        # 返り値が２つあるので、割り当てる場合には左辺の変数を ',' で区切る \n",
    "        b0hat, b1hat = my_sim_manual(n,b0=b0,b1=b1,su=su)\n",
    "\n",
    "        b0hat_arr[i] = b0hat  # i番目に定数項の推定値をarrに追加\n",
    "        b1hat_arr[i] = b1hat  # i番目にスロープ係数の推定値をarrに追加\n",
    "        \n",
    "    return b0hat_arr, b1hat_arr  # N回ループ終了後にarrayをタプルとして返すオブジェクト"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "標本の大きさを`50`，シミュレーションの回数を`100_000`としてシミュレーションを始める。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "b0hat, b1hat = my_sim_loop(n=50,N=100_000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "```{note}\n",
    "初めて`my_sim_loop()`を実行すると`Numba`はコンパイルと呼ばれる高速化準備プロセスを経るため，`@njit`がない場合と比べて実行スピードは早くない。２回目の実行から`Numba`の効果がフルに発揮されることになる。\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "````{hint}\n",
    "`Python`の`for`ループはコンパイル言語（例えば，`C`）と比べると実行速度が遅い。それを少しでも克服しようというのが`Numba`である。一方，`for`ループ自体を使わずに`Numpy`の特徴であるベクトル演算を使うと格段に実行速度が速くなる。例えば，`my_sim_loop`の代わりに次の関数を試してみよう。\n",
    "```\n",
    "def my_sim_vectorize(n,N,b0=1.0,b1=1.0,su=1.0):\n",
    "    \n",
    "    # 一度に全てのランダム変数を生成し，(n,N)のarrayに変換する\n",
    "    x = np.random.normal(loc=4, scale=1, size=(n,N))\n",
    "    u = np.random.normal(loc=0, scale=su, size=(n,N))\n",
    "    y = b0 + b1*x + u\n",
    "\n",
    "    # 平均の計算\n",
    "    x_mean = x.mean(axis=0)\n",
    "    y_mean = y.mean(axis=0)\n",
    "\n",
    "    # 分散の計算\n",
    "    x_var = x.var(axis=0)\n",
    "    y_var = y.var(axis=0)\n",
    "\n",
    "    # 平均からの乖離を計算\n",
    "    x_mean_dev = x - x_mean\n",
    "    y_mean_dev = y - y_mean\n",
    "\n",
    "    # 共分散の計算\n",
    "    yx_dev = y_mean_dev * x_mean_dev\n",
    "    yx_cov = yx_dev.mean(axis=0)\n",
    "    \n",
    "    # 推定値の計算\n",
    "    b1 = yx_cov / x_var\n",
    "    b0 = y_mean - b1*x_mean\n",
    "    \n",
    "    return b0, b1\n",
    "```\n",
    "この関数には`for`ループはなく，戻り値`b0`と`b1`は`N`個の要素からなる推定値のarrayとなる。実行速度を比べるために次の２つのコードを別々のコードセルで実行してみよう。\n",
    "```\n",
    "%timeit my_sim_loop(n=50,N=1_000_000)\n",
    "```\n",
    "```\n",
    "%timeit my_sim_vectorize(n=50,N=1_000_000)\n",
    "```\n",
    "ここで`%timeit`は実行速度を測るマジック・コマンドと呼ばれるコードである。`my_sim_vectorize`が２倍近く速くなることが確認できるだろう。`Numba`のようなパッケージを導入しなくても，ベクトル演算を使う方が計算の時間を短縮できることを覚えておこう。\n",
    "````"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "$\\hat{\\beta}_0$の分布を図示する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "plt.hist(b0hat,bins=60, edgecolor='white')\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "$\\hat{\\beta}_1$の分布を図示する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "plt.hist(b1hat,bins=60, edgecolor='white')\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`y`の理論値を図示するが，図示の時間短縮のために`N`回のシミュレーション中の最初の`r`回の回帰線だけを表示する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "r = 1000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`r`回の回帰線を`for`ループを使って重ねて図示してみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "xx = np.linspace(2,6,50)\n",
    "\n",
    "for i in range(r):  # r回のループで図示\n",
    "    yhat = b0hat[i] + b1hat[i]*xx  # yの予測値\n",
    "    plt.plot(xx, yhat, 'grey')  # yの予測値を図示\n",
    "    \n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.ylim(0,10) # 縦軸の表示範囲を設定\n",
    "plt.title(f'{r} rounds of OLS simulation')\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "```{admonition} コードの説明\n",
    "上のコード・ブロックの最後から２行目に`f'{r} rounds of OLS simulation'`とある。文字列の前に`f`を置き、文字列の中で上で定義した変数を`{}`で囲むとその変数の値を代入することができる。これを`f-string`と呼ぶ。\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "**＜シミュレーションの含意＞**\n",
    "\n",
    "まず上の２つのヒストグラムをみると，値にばらつきがある。これは推定量（$\\hat{\\beta}_0$と$\\hat{\\beta}_1$）が確率変数であるためである。その実現値は，真の値（$\\beta_0=1.0$と$\\beta_1=1.0$）を中心に左右対象に分布している。即ち，推定量は不偏性を満たしている。次の章で簡単に説明するが，上で行ったシミュレーションは不偏推定量に必要な仮定を満たしているためである。\n",
    "\n",
    "これを踏まえて，回帰直線の図を考えよう。同じ回帰式を使い`N`回のシミュレーションを行った後，最初の`r`個の回帰直線を重ねて描いているためパイプのような形になっている。一方で，実際の実証分析ではこのような繰り返しの試行は不可能であり多くの場合は一度だけである。言い換えると，実証研究では`N`個の回帰直線の中から１つだけを結果として「ランダムに選んでいる」のと等しい。この意味を更に考えるために，（標本の大きさは所与として）たまたま２組みのデータセットが手に入り，別々の回帰分析を行ったとしよう。その場合，パイプの幅が広ければ広いほど，２つの回帰直線の位置は大きく異なる可能性が高くなり，正確性が低い結果になる。逆に，パイプの幅が狭ければ狭いほど，２つの結果は似たものになり，より正確な推定結果となる。このように実証分析ではパイプの幅が可能な限り小さくなる工夫が必要になる。その１つが標本の大きさ`n`である。`n`が増えるとパイプは細くなり，極限$n\\rightarrow\\infty$（即ち，標本の大きさと母集団の大きさの差が限りなくゼロに近づく）ではパイプの形が直線（母集団回帰式）に近づくことになる。このように，実証分析において標本の大きさは非常に重要な役割を果たすことになる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "```{warning}\n",
    "* このシミュレーションの場合，推定量の**一致性**が満たされているため，パイプが母集団回帰式に収束することになる。しかし一般的に必ずそうなるとは限らない。ある仮定のもとでのみ推定量は一致性を満たすことになる。\n",
    "* 推定方法によっては，不偏性は満たさないが一致性を満たす推定量も存在する。また，エギゾチックな例になるがその逆もあり得る。\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "----- Give it a try! -----\n",
    "\n",
    "標本の大きさ`n`を変えて試してみよう。例えば，10, 30, 100, 500, 1000。"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernel_info": {
   "name": "py4etrics"
  },
  "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.14.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": 4
}