9_Multiple_Regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(chap:9-multiple)=\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 seaborn as sns\n",
    "import statsmodels.formula.api as smf\n",
    "import py4macro\n",
    "import wooldridge\n",
    "\n",
    "from numba import njit\n",
    "from pandas.plotting import scatter_matrix\n",
    "from scipy.stats import norm, uniform, gaussian_kde, multivariate_normal\n",
    "from statsmodels.stats.outliers_influence import variance_inflation_factor as vif\n",
    "\n",
    "# 警告メッセージを非表示\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## 説明"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "説明変数が複数の重回帰分析（Multiple Regression）を考える。\n",
    "\n",
    "$$y_i=\\beta_0+\\sum_{j=1}^k\\beta_j x_{ij} + u_i$$\n",
    "* $i=1,2,...n$：観測値のインデックス\n",
    "* $j=1,2,...k$：説明変数の数\n",
    "\n",
    "（注意）\n",
    "\n",
    "$x_1$と書く場合，$x_{i1},x_{i1},...x_{n1}$を表す（第１番目の説明変数）。\n",
    "\n",
    "一般的に$x_j$と書く場合，$x_{ij},x_{ij},...x_{nj}$を表す（第$j$番目の説明変数）。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "## `statsmodels`を使う"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "### 推定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "単回帰分析で使った同じデータセット`wage1`を使い，賃金と教育の関係を再考する。[前章](chap:8-simple)と異なる点は，複数の説明変数を導入し重回帰分析をおこなうことである。追加的な説明変数として次の変数を加える。\n",
    "* `exper`：潜在的な経験年数（＝年齢ー最高学歴までの年数ー６）\n",
    "* `tenure`：現在の職場での在職年数\n",
    "\n",
    "`statsmodels`を使う上で単回帰と異なる点は，それぞれの説明変数の前に`+`を追加するだけである。\n",
    "\n",
    "まず必要な変数だけを取り出し変数`wage1`に割り当てる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "wage1 = wooldridge.data('wage1').loc[:,['wage', 'educ', 'tenure', 'exper']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true,
    "tags": [
     "output_scroll"
    ]
   },
   "outputs": [],
   "source": [
    "wooldridge.data('wage1',description=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "被説明変数に対数賃金を使うために、新たな列`wage_log`を追加する"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "wage1['wage_log'] = np.log(wage1['wage'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "回帰式では追加の説明変数の前に`+`を付ける。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "formula_1 = 'wage_log ~ educ + tenure + exper'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "推計結果を変数`res_1`に割り当てる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "res_1 = smf.ols(formula_1, data=wage1).fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "[単回帰分析](sec:8-statsmodels)で説明したが，次のコードで結果を表にまとめて表示することができる。（関数`print()`を使っているが使わなくても同じ情報が表示される。）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "print(res_1.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "### 属性とメソッド"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "推定結果`res_1`の属性とメソッドは`dir()`や`py4macro.see()`で確認できる。以下では`py4macro.see()`を使う。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "output_scroll"
    ]
   },
   "outputs": [],
   "source": [
    "py4macro.see(res_1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "単回帰分析と同様に、係数の推定値は属性`.params`で取得できる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_1.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`res_1`は`Series`なので、インデックスを使って`educ`の係数だけを取り出す場合は次のようにする。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_1.params.iloc[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "次に，メソッド`predict()`を使って予測値を計算してみよう。次の例では，`educ`, `tenure`, `exper`の平均値での`np.log(wage)`を計算する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wm = wage1.mean()  # 変数の平均の計算\n",
    "wm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`wm`は`Series`なので、変数の行ラベルを使い`wm['educ']`で`educ`の平均値だけを取り出すことができる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wm['educ']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`predict()`に渡す引数としては，説明変数の辞書を作るのが一番簡単であろう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z = {'educ':wm['educ'],\n",
    "     'tenure':wm['tenure'],\n",
    "     'exper':wm['exper']}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_1.predict(z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "３つの変数が平均値を取る場合の平均時給の対数値は約1.62である。\n",
    "\n",
    "複数の予測値を計算するには、辞書の値に`list`や`array`を書き込めば良い。次の例では、`tenure`と`exper`を平均値に固定し、`educ`を平均値、平均値-1、平均値+1に変化させて予測値を計算している。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "educ_mean = wm['educ']\n",
    "tenure_mean = wm['tenure']\n",
    "exper_mean = wm['exper']\n",
    "\n",
    "z2 = {'educ':[educ_mean-1, educ_mean, educ_mean+1],\n",
    "                  'tenure':[tenure_mean]*3,\n",
    "                  'exper':[exper_mean]*3}\n",
    "\n",
    "educ_return = res_1.predict(z2)\n",
    "educ_return"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "３つの変数が平均の値を取る場合，平均賃金の対数値は1.623268であり，上の計算と同じになることが確認できる。教育年数が1年減少すると予測値は1.531239であり，教育年数が1年増えると予測値は1.715297となる。この結果を使い，教育年数が一年変化する場合の影響を計算してみよう。まず教育年数が1年減少すると場合を計算しよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "educ_return[0]-educ_return[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "次に教育年数が1年増える場合を計算しよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "educ_return[2]-educ_return[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "これらの数字の絶対値は同じである。この結果は偶然ではない。他の要因を固定し教育年数が1年変化（減少もしくは増加）した場合，平均時給の対数値がどれだけ変化したかを計算している。これは正しく教育の収益率である。この点を数式で確認するために，変数$x$の値が小さい場合（例えば，0.02），次の近似が成立することを思い出そう。\n",
    "\n",
    "$$\n",
    "\\log(x+1)\\approx x\n",
    "$$ (eq:9-approx)\n",
    "\n",
    "また教育年数が$h$の場合の平均時給を$w(h)$としよう。教育年数が一年増えた場合の教育の収益率は$\\dfrac{w(h+1)-w(h)}{w(h)}$となり，[](eq:9-approx)を使うと次のように書くことができる。\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\frac{w(h+1)-w(h)}{w(h)}\n",
    "&\\approx\\log\\left(\\frac{w(h+1)-w(h)}{w(h)}+1\\right)\\\\\n",
    "&=\\log\\left(\\frac{w(h+1)}{w(h)}\\right)\\\\\n",
    "&=\\log(w(h+1))-\\log(w(h))\n",
    "\\end{align*}\n",
    "$$ (eq:9-approx2)\n",
    "\n",
    "最後の等号が上の計算で使った式である。$h+1$を$h-1$に変えると教育年数が1年減少した場合を考えることができる。\n",
    "\n",
    "もう一点付け加えるために，上の回帰分析の変数`educ`の係数の解釈を思い出そう。教育年数が一年増えた場合，平均賃金の対数値はどれだけ変化するかを示しており，教育の収益率を表している。従って，ここで計算した値と同じになるはずである。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_1.params['educ']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "（コメント）\n",
    "\n",
    "`predict()`の引数を省略すると予測値を返す属性`.fittedvalues`と同じ値を返す。それを確かめるために次のコードを評価してみる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "( res_1.predict() == res_1.fittedvalues ).all()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{admonition} コードの説明\n",
    "* `res_1.predict() == res_1.fittedvalues`の意味\n",
    "    * `res_1.predict()`と`res_1.fittedvalues`はSeriesを返すので，`==`で比較することにより，それぞれの対応する要素の値が等しければ`True`，異なれば`False`になるSeriesを返す。\n",
    "* `( res_1.predict() == res_1.fittedvalues )`の`()`は、その中のコードを先に実行することを指定する（数学の`()`と同じ）。\n",
    "* `all()`は返されたSeriesの値が全て`True`である場合は`True`を、１つでも`False`があれば`False`を返す。\n",
    "```\n",
    "\n",
    "`True`が返されているので、値は全て同じだと確認できた。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "## 変数の変換"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "上の重回帰分析では，回帰式の`wage`を対数化した変数を新たに作成し計算を行なっているが、回帰式の中で変数の変換を直接指定することもできる。以下で説明する方法は被説明変数・説明変数の両方に使用可能である。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true,
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "### `Numpy`の関数を使う方法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`Numpy`の対数関数を使って，回帰式の中で直接書き換えることができる。`wage`を対数化する例を考える。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "formula_2 = 'np.log(wage) ~ educ + tenure + exper'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "res_2 = smf.ols(formula_2, data=wage1).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "res_2.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true,
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "### `I()`を使う方法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "次に`exper`の二乗を回帰式に入れるケースを考えよう。この場合、`np.square(exper)`でも良いが、`I()`の`()`の中に直接式を書くことが可能となる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "formula_3 = 'np.log(wage) ~ educ + tenure + exper + I(exper**2)'\n",
    "res_3 = smf.ols(formula_3, data=wage1).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "res_3.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "この方法であれば、`Numpy`に無い関数も直接書くことができる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true,
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "### `def`を使う方法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "より複雑な関数であれば`def`を使い、関数を定義し、それを回帰式の中で使う方が良いかも知れない。`exper`の二乗を回帰式に入れるケースを考えよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def myfunc(x):\n",
    "    return x**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "formula_4 = 'np.log(wage) ~ educ + exper + tenure + myfunc(exper)'\n",
    "res_4 = smf.ols(formula_4, data=wage1).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "res_4.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true,
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "### `Q()`について"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`Q()`は変数の変換に使うわけではないが、その用途を覚えておくと役に立つだろう。\n",
    "\n",
    "`Python`には予約語（reserved words）と呼ばれるものがり、変数名には使えない。以下がその例である。\n",
    "* `def`と`return`：関数を定義する場合に使う\n",
    "* `for`と`while`：ループに使う\n",
    "* `import`と`from`：パッケージのインポート\n",
    "* `if`と`elif`と`else`：`if`文に使う\n",
    "\n",
    "また`DataFrame`の列ラベルに予約語を文字列として使うことも避けるべきであろう。しかし、どうしても使う場合もあるかもしれない（例えば、他のプログラム用に作成されてデータを読み込む場合）。だが、予約語をそのまま回帰式に書くとエラーが発生する。その際に役に立つのが`Q()`である。例えば、`df`という`DataFrame`に「利子率」の列があり，列ラベルとして`return`が設定されているとしよう。`Q()`の引数として`return`を文字列で書くことにより、通常どおりのコードを書くことができる。\n",
    "\n",
    "```\n",
    "formula = 'Q(\"return\") ~ X\n",
    "res = ols(formula, data=df).fit()\n",
    "```\n",
    "\n",
    "また`123go`のように数字から始まる列ラベルもあるかもしれない。この場合も回帰式でそのまま使うとエラーが発生するが、`Q()`の引数として文字列で指定すればが使えるようになる。\n",
    "```\n",
    "formula = 'Q(\"return\") ~ X + Q(\"123go\")'\n",
    "res = ols(formula, data=df).fit()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "tags": [
     "remove-cell"
    ]
   },
   "source": [
    "## 係数の解釈"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true,
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "(sec:9-interpretation)=\n",
    "## 係数の解釈"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "次の推定式を考えよう。\n",
    "\n",
    "$$y=\\beta_0+\\beta_1x_1+\\beta_2x_2+\\beta_3x_3+u$$\n",
    "\n",
    "＜解釈１：Ceteris Paribus＞\n",
    "\n",
    "$x_1$以外の変数を固定して$x_1$だけを少しだけ変化させる。これにより，$y$に対する$x_1$の影響だけを考えることが可能となり、次の記号で表す。\n",
    "\n",
    "$$\\dfrac{\\Delta y}{\\Delta x_1}=\\beta_1$$\n",
    "\n",
    "ここで\n",
    "* $\\Delta x_1=x_1$の変化（１単位の変化と解釈する） \n",
    "* $\\Delta y=y$の変化 \n",
    "\n",
    "$\\beta_1$は、$x_1$が一単位増加した場合の$y$の変化を示している。重要な点は，$x_1$以外の変数は固定されていることである。これをCeteris Paribus (ラテン語; \"other things being equal\"という意味)と呼ぶ。\n",
    "\n",
    "＜解釈２：影響の除去＞\n",
    "\n",
    "まず$y$と$x_1$の関係を捉える$\\hat{\\beta}_1$を考えよう。もし推定式に$x_2$と$x_3$が入っていなければ，$\\hat{\\beta}_1$は$x_2$と$x_3$の影響を拾うことになる（$x_1$が$x_2$と$x_3$のそれぞれと完全に無相関でない場合）。即ち，$x_2$と$x_3$が入ることにより，$\\hat{\\beta}_1$はこれらの変数の影響を「除去」した後に残る$y$と$x_1$の関係を捉えていると考えることができる。同じように，$\\hat{\\beta}_j$は$j$以外の説明変数の影響を取り除いた後に残る$y$と$x_j$の関係を示している。\n",
    "\n",
    "---\n",
    "`wage1`を使い以下を推定しよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "formula_1a = 'wage ~ educ + tenure + exper'\n",
    "res_1a = smf.ols(formula_1a, data=wage1).fit()\n",
    "res_1a.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`tenure`と`exper`を固定して，`educ`を一単位増加させるとしよう。その効果は\n",
    "\n",
    "$$\\dfrac{\\Delta\\text{wage}}{\\Delta\\text{educ}}\\approx 0.599$$\n",
    "\n",
    "即ち，教育期間が一単位（ここでは一年）増えると賃金は0.599ドル増加することを示している。\n",
    "\n",
    "---\n",
    "次に`formula_1`の回帰式を使って推定した`res_1`を考えよう。この場合，\n",
    "\n",
    "$$\\dfrac{\\Delta\\ln(\\text{wage})}{\\Delta\\text{educ}}\\approx 0.092$$\n",
    "\n",
    "対数の性質を使い左辺は次のように書き換えることができる。\n",
    "\n",
    "$$\\dfrac{\\Delta\\ln(\\text{wage})}{\\Delta\\text{educ}}\n",
    "=\\dfrac{\\Delta\\text{wage}/\\text{wage}}{\\Delta\\text{educ}}\n",
    "=\\dfrac{\\%\\Delta\\text{wage}}{\\Delta\\text{educ}}\n",
    "$$\n",
    "\n",
    "ここで，$\\%\\Delta\\text{wage}=$％で表した\n",
    "`wage`の変化を意味する（[](eq:9-approx2)を参照）。従って，教育期間が一単位（ここでは一年）増えると賃金は約9.2%増加することを示している。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## ガウス・マルコフ定理"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### 母集団回帰式と標本回帰式"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "今までは回帰式に基づいて`Python`を使い係数の推定方法を説明したが，ここではOLS推定量の望ましい性質とはなにか，望ましい性質の条件はなにかを考える。そのために，まず母集団回帰式と標本回帰式について説明する。\n",
    "\n",
    "計量経済学では，データを使い**経済理論（経済学的な考え）**\n",
    "に基づいた仮説を検証する。例として「消費は家計の所得$x_1$が高くなると増え，子どもの数$x_2$とも正比例する」という仮説を考えよう。所得と子どもの数に依存する消費の部分だけを$y^*$と書くと，一般型の消費関数は$y^*=f\\left(x_1,x_2\\right)$となる。更に，線形を仮定すると次式となる。\n",
    "\n",
    "$$y^*=\\beta_0+\\beta_1x_1+\\beta_2x_2\\qquad{(式１)}$$\n",
    "\n",
    "例えば，所得30万円，子ども2人の家計の場合、消費は$\\beta_0+\\beta_1\\cdot\\left(\\text{30万円}\\right)+\\beta_2\\cdot\\left(\\text{2人}\\right)$となる。所得 $x_1$が変化すると消費$y^*$も変化することが分かる。このことは$\\beta_2=0$の単回帰の場合（消費が子供の数に依存しない場合）を考えればもっと分かりやすいだろう。その場合，(式１)は\n",
    "\n",
    "$$y^*=\\beta_0+\\beta_1x_1$$\n",
    "\n",
    "となり，横軸に$x_1$，縦軸に$y^*$を置く図では右上がりの直線となる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "(式１)は**平均**の家計の行動を表していると考えることができる。従って，同じ所得と子どもの数の家計でも好み等の要因が違うと消費に違いがでる。その違いを（平均ゼロの）誤差項$u$として捉えると，全ての家計の消費 $y$は次式で表すことができる。\n",
    "\n",
    "$$y=y^*+u\\qquad{(式２)}$$\n",
    "\n",
    "* $y^*$：平均の家計の消費であり，$x_1$と$x_2$のみに依存する部分\n",
    "* $u$：平均の家計と異なる家計の消費であり，平均家計との違いを捉える部分\n",
    "\n",
    "(式１)を(式２)に代入すると\n",
    "\n",
    "$$y=\\beta_0+\\beta_1x_1+\\beta_2x_2+u\\qquad{(式３)}$$\n",
    "\n",
    "となる。この式は仮説に基づいた関係式であり，**母集団回帰式**と呼ばれる。また(式３)を使い$x_1$と$x_2$を所与として期待値を計算すると次式を得る。\n",
    "\n",
    "$$\\text{E}\\left(y|x_1,x_2\\right)=y^*\\qquad\\because\\;\\text{E}(u)=0$$\n",
    "\n",
    "これは「$x_1$と$x_2$を所与とする」という条件の下での期待値なので**条件付き期待値**（conditional expectations）と呼ばれる。この式は，$y$の条件付き期待値は母集団回帰式の$x_1$と$x_2$に依存する部分である$y^*$と等しいことを示している。即ち，（式１）は母集団で平均として成立する**母集団回帰線**を表している。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "```{note}\n",
    "* 母集団回帰式(式３)は理論的な仮定に基づいて導出したため，それが正しいかどうかは事前に分からない。従って，経済理論に基づき母集団回帰式の妥当性を検討する必要がある。\n",
    "* 母集団パラメータ$\\beta_j$, $j=0,1,2$と誤差項$u$は**観測できない**。\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "母集団回帰式(式３)が成立するという**仮定**のもとで母集団パラメータ$\\beta_j$，$j=0,1,2$を推定することになる。そのためには，\n",
    "* $\\left\\{y_i,x_{i1},x_{i2}\\right\\}$, $i=1,2,..,n$ の標本を収集する。\n",
    "* 次式に重回帰分析の手法を適用しOLS推定値を計算する。\n",
    "\n",
    "    $$y_i=\\beta_0+\\beta_1x_{i1}+\\beta_2x_{i2}+u_i,\\qquad i=1,2,..,n\\qquad{(式４)}$$\n",
    "\n",
    "母集団回帰式(式３)と区別するために，この式には添え字$i$があり**標本回帰式**と呼ぶ。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "```{note}\n",
    "通常，標本の大きさは母集団のそれよりも小さい。従って，パラメータのOLS推定値$\\hat{\\beta}_j$には誤差が発生する。即ち，必ずしも$\\hat{\\beta}_j=\\beta_j$, $j=0,1,2$とはならない。\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "一方で，上で述べたようにOLS推定量には何らかの「望ましい性質」がある。\n",
    "1. 「望ましい性質」の前提となる条件は何なのか？ \n",
    "1. OLS推定量の「望ましい性質」とは何なのか？\n",
    "\n",
    "以下では１から考える。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true,
    "tags": [
     "remove-cell"
    ]
   },
   "source": [
    "### 仮定と望ましい性質"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "(sec:9-assumptions)=\n",
    "### 仮定と望ましい性質"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "* 仮定１：Linear in Parameters（母集団回帰式の線形性）\n",
    "    * 母集団のモデルは次式のようにパラメータに関して線形である。\n",
    "    \n",
    "        $$\n",
    "        y=\\beta_0+\\beta_1x_1+\\beta_2x_2+\\cdots+\\beta_kx_k+u\n",
    "        $$\n",
    "\n",
    "    * $y=\\beta_0+\\beta_1x_1+\\beta_2x_1^2$は線形性を満たす。\n",
    "    * $y=\\beta_0+\\dfrac{\\beta_1}{1+\\beta_2x_1}+\\beta_3x_2$は線形性を満たさない。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "* 仮定２：Random Sampling（標本の無作為抽出）\n",
    "    * 無作為に$\\left\\{y_i,x_{i1},,x_{i2},...,,x_{ik}\\right\\}$, $i=1,2,...,n$が抽出される。\n",
    "    * 多くの横断面データでは成立すると考えられている（時系列データでは満たされない場合が多い）。\n",
    "    * この仮定により\n",
    "    \n",
    "        $$\\text{Corr}\\left(u_i,u_j|X\\right)=0,\\;i\\neq j,\\qquad X=\\{x_1,x_2,...x_k\\}$$\n",
    "        \n",
    "        となる。つまり，説明変数$X$を所与として$u_i$と$u_j$は何の関係もないという意味である。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "* 仮定３：No Perfect Collinearity（母集団回帰式の説明変数の非完全多重共線性）\n",
    "    * 説明変数$x_j$が他の説明変数の完全な線形結合（例えば，$x_2=1+2x_1$や$x_1=x_2+5x_4$）として表すことができないという意味である。もしこの仮定が満たされなければ$x_j$は何の追加的な情報を持たなくなる。\n",
    "    * 定数項以外の説明変数が全く変化しないことも排除している（例えば，$x_j=3$）。このような変数は定数項の単なるスケールアップ・バージョンであり，追加的な情報はない。（定数項以外の説明変数は変化してこそ$y$を説明できる。）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "* 仮定４：Zero Conditional Mean（母集団回帰式の誤差項の条件付き期待値は`0`）\n",
    "\n",
    "    $$\\text{E}\\left(u|X\\right)=0,\\qquad X=\\{x_1,x_2,...x_k\\}$$\n",
    "    \n",
    "    * $X$はランダム変数であるが，$X$を所与とすると$u$の平均はゼロになるという仮定である。即ち，ランダム変数である$X$がどのような値を取ったとしても，$u$の平均は0だということである。$X$は$u$に対して何の情報も影響力も全くないため，$X$と$u$が線形・非線形な関係が無いことを意味する。即ち，$X$と$u$は独立性を満たす。\n",
    "    * $X$はランダム変数であるため，$X$が変化する場合の$\\text{E}\\left(u|X\\right)$の平均を考えることができる。即ち，$\\text{E}\\left[\\text{E}\\left(u|X\\right)\\right]=0$である。外側の$E$は$X$に対する期待値演算子であり，$E$が入れ子になっているのでiterated expectationsと呼ばれる。$X$がどのような値を取ったとしても仮定４のもとでは$0$になるので等号が成立することになる。更に，$X$に対しての期待値が付け加えられているので，$X$を省いて$\\text{E}\\left[\\text{E}\\left(u|X\\right)\\right]=\\text{E}(u)=0$となる。$X$の条件をつけなくとも$u$の平均はゼロになることを示している。\n",
    "    * $X$と$u$には線形の依存関係がないので$\\text{Cov}(X,u)=0$が成立する。しかし，$\\text{Cov}(X,u)=0$（無相関）は$\\text{E}\\left(u|X\\right)=0$（独立性）を意味しないことを覚えておこう。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "仮定１〜４の下で母集団パラメータのOLS推定量は不偏性（unbiasedness）を満たす。\n",
    "\n",
    "$$\\text{E}\\left(\\hat{\\beta}_j\\right)=\\beta_j,\\qquad j=0,1,2,...,k$$\n",
    "\n",
    "（解釈）\n",
    "\n",
    "標本の大きさが$n$の標本を使い回帰分析をおこない$\\hat{\\beta}_j$を得たとしよう。さらに同じ標本の大きさ$n$の異なる標本を使い回帰分析を合計$N$回おこない$\\hat{\\beta}_j$を計算したとしよう。不偏性とは$N$が十分に大きい場合（$N\\rightarrow\\infty$），$N$個ある$\\hat{\\beta}_j$の平均は$\\beta_j$と等しいという性質である。\n",
    "\n",
    "（注意点）\n",
    "\n",
    "この結果は標本の大きさ$n$に依存しない。即ち，大標本（$n\\rightarrow\\infty$）である必要はない。むしろ小標本（$n<\\infty$）で重要視される小標本特性（データが少ない場合の特性）と考えるべきである。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "* 仮定５：Homoskedasticity（均一分散; 母集団回帰式の説明変数は誤差項の分散に影響を与えない）\n",
    "\n",
    "$$\\text{Var}\\left(u|X\\right)=\\sigma^2$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "仮定１〜５をガウス・マルコフ仮定（横断回帰分析のGM仮定）と呼び，この仮定の下では以下が成立する。\n",
    "1. 誤差分散の推定量 $\\hat{\\sigma}^2$の平均は母集団の値と等しい（誤差分散の不偏推定量）。\n",
    "\n",
    "    $$\n",
    "    \\text{E}\\left(\\hat{\\sigma}^2\\right)=\\sigma^2\n",
    "    $$\n",
    "\n",
    "    ここで，$\\hat{\\sigma}^2=\\dfrac{\\sum_{i=1}^n\\hat{u}_i^2}{n-k-1}=\\dfrac{SSR}{n-k-1}$\n",
    "    であり，$\\hat{u}_i=y_i-\\hat{y}_i$はOLS残差であり，$SSR$は残差の二乗平方和。$\\hat{\\sigma}$にはさまざまな呼称があり，the Standard Error of Regression, the Standard Error of the Estimate, the Root Mean Squared Errorがある。\n",
    "\n",
    "2. OLS係数の（標本）分散の推定量（$\\sigma^2$に$\\hat{\\sigma}^2$が代入されている）\n",
    "\n",
    "    $$\n",
    "    \\widehat{\\text{Var}}\\left(\\hat{\\beta}_j\\right)\n",
    "    =\\dfrac{\\hat{\\sigma}^2}{SST_j\\left(1-R_j^2\\right)}\n",
    "    =\\dfrac{1}{n-1}\\cdot\n",
    "    \\dfrac{\\hat{\\sigma}^2}{\\text{Var}(x_j)}\\cdot\n",
    "    \\dfrac{1}{1-R_j^2}\n",
    "    $$\n",
    "\n",
    "    ここで，$\\text{Var}(x_j)=(n-1)^{-1}\\sum_{i=1}^n\\left(x_{ij}-\\overline{x}_j\\right)$であり，$R_j^2$は$x_j$を他の全ての説明変数に回帰分析する際の決定係数である。\n",
    "    $\\widehat{\\text{Var}}\\left(\\hat{\\beta}_j\\right)$はOLS推定量$\\hat{\\beta}_j$の精度を捉え，低ければより高い精度を意味する。（「なぜ分散？」と思う場合は前章のシミュレーション１で，標本のランダム抽出を行うたびに，$\\hat{\\beta}$の値が変わったことを思い出そう。）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "パラメータの分散の推定量の式から次のことが分かる。\n",
    "* $\\dfrac{1}{n-1}$: 標本の大きさが大きくなると$\\hat{\\beta}_j$の精度を高める。\n",
    "* $\\dfrac{\\hat{\\sigma}^2}{\\text{Var}(x_j)}$: 残差の分散に対して説明変数$x_j$の分散が大きければ，$\\hat{\\beta}_j$の精度を高める。\n",
    "* $R_j^2$は多重共線性が高くなると大きくなる。従って，多重共線性が高いと$\\hat{\\beta}_j$の精度を低下させる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "$\\widehat{\\text{Var}}\\left(\\hat{\\beta}_j\\right)$のルートは**OLS推定量の標準誤差**と呼ぶ（パラメータの推定値の標準偏差）。\n",
    "\n",
    "$$\n",
    "\\text{se}\\left(\\hat{\\beta}_j\\right)\n",
    "=\\sqrt{\\widehat{\\text{Var}}\\left(\\hat{\\beta}_j\\right)}\n",
    "=\\dfrac{1}{\\sqrt{n-1}}\\cdot\n",
    "\\dfrac{\\hat{\\sigma}}{\\text{sd}(x_j)}\\cdot\n",
    "\\dfrac{1}{\\sqrt{1-R_j^2}}\n",
    "$$\n",
    "\n",
    "この結果は推定結果の検定や信頼区間を計算する際に必要となる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "**ガウス・マルコフ定理**（望ましい性質）\n",
    "\n",
    "GM仮定１〜５のもとで`OLS`推定量はB.L.U.E（Best Linear Unbiased Estimators; 最良不偏推定量）\n",
    "* 不偏性（unbiasedness）：\n",
    "$\\text{E}\\left(\\hat{\\beta}_j\\right)=\\beta_j$\n",
    "* 効率性（efficiency）： $\\text{Var}\\left(\\hat{\\beta}_j\\right)$は線形推定値の中で最小"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "---\n",
    "上の回帰分析結果`res_1`を使って，ここで定義した推定値の取得方法を説明する。\n",
    "\n",
    "**誤差分散$\\hat{\\sigma}^2$の推定値**\n",
    "\n",
    "`res_1`の属性`.scale`で取得可能。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "res_1.scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "SSR = res_1.ssr  # 残差の二乗平方和\n",
    "n = res_1.nobs\n",
    "k = res_1.df_model\n",
    "SSR/(n-k-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "```{note}\n",
    "上の計算に属性`.ssr`が出てくる。インターネットで`statsmodels result`検索してみよう。[Google](https://www.google.com/search?client=safari&rls=en&q=statsmodels+result&ie=UTF-8&oe=UTF-8)の検索結果の[一番上のサイト](https://www.statsmodels.org/devel/generated/statsmodels.regression.linear_model.RegressionResults.html)は該当する`statsmodels`マニュアルとなっており、ウェブページの下の方に行くと`.ssr`を見つけることができるはずだ。\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "**パラメータの標準誤差 $\\text{se}\\left(\\hat{\\beta}_j\\right)$**\n",
    "\n",
    "`res_1`のメソッド`.bse`で取得可能。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "res_1.bse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "パラメータの標準誤差は$se$（standard error）だが属性は`bse`となっている。これはパラメータを$b$で表しており，$b$の$se$として覚えれば良いだろう。"
   ]
  },
  {
   "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": [
    "`statsmodels`は非常に優れたパッケージであり，複雑な計算を簡単なコードで実行してくれる。しかしここでは`Numpy`の関数を使い，重回帰分析における推定値を「手計算」で求める。目的は２つある。第１に，`statsmodels`の裏でどのような計算が行われているかを理解する。第２に，今後シミュレーションをおこなうが，「手計算」のコードを使うと`Numba`を使い計算速度を格段とアップすることが可能となる。\n",
    "\n",
    "次の回帰式を考えよう。\n",
    "\n",
    "$$\n",
    "y_i=\\beta_0+\\beta_1x_{1i}+\\beta_2x_{2i}+u_i,\\quad\n",
    "i=1,2,...,n\n",
    "$$\n",
    "\n",
    "これを行列式で表すと次のようになる。\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\begin{bmatrix}\n",
    "y_1 \\\\\n",
    "y_2 \\\\\n",
    "\\vdots \\\\\n",
    "y_n\n",
    "\\end{bmatrix}\n",
    "=&\n",
    "\\begin{bmatrix}\n",
    "1,x_{11},x_{21} \\\\\n",
    "1,x_{12},x_{22} \\\\\n",
    "\\vdots\\vdots\\vdots \\\\\n",
    "1,x_{1n},x_{2n}\n",
    "\\end{bmatrix}\n",
    "\\cdot\n",
    "\\begin{bmatrix}\n",
    "\\beta_0 \\\\\n",
    "\\beta_1 \\\\\n",
    "\\beta_2\n",
    "\\end{bmatrix}\n",
    "+\n",
    "\\begin{bmatrix}\n",
    "u_1 \\\\\n",
    "u_2 \\\\\n",
    "\\vdots \\\\\n",
    "u_n\n",
    "\\end{bmatrix}\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "更に簡単に書くと\n",
    "\n",
    "$$Y=XB+U$$\n",
    "\n",
    "となる。これを使うと$B$の推定量は\n",
    "\n",
    "$$\n",
    "\\hat{B}=\\left(X^T\\cdot X\\right)^{-1}\\cdot X^T\\cdot Y\\qquad\\qquad\\text{(*)} \n",
    "$$\n",
    "\n",
    "で与えられる。以下ではシミュレーションを使い式（\\*）の計算に使うコードを説明する。\n",
    "\n",
    "まず母集団回帰式の係数と標本の大きさを設定する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "b0 = 1  # 定数項\n",
    "b1 = 2  # x1の係数\n",
    "b2 = 3  # x2の係数\n",
    "n = 30  # 標本の大きさ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "母集団からの標本を生成する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "x1 = np.random.normal(loc=4.0, scale=2.0, size=n)  # (1)の説明\n",
    "\n",
    "x2 = np.random.uniform(low=1.0, high= 10.0, size=n)  # (2)の説明\n",
    "\n",
    "u = np.random.normal(loc=0, scale=1.0, size=n)  # (1)の説明\n",
    "\n",
    "y = b0 + b1*x1 + b2*x2 + u\n",
    "\n",
    "c = np.ones(n)  # (3)の説明"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "1. `Numpy`には乱数用のモジュール`random`があり，`normal`は標準正規分布からの乱数を発生させる関数である。引数の設定方法は`scipy.stats`の`norm.rvs()`と同じである。`x1`は平均`4.0`と標準偏差`2.0`の正規分布に従うことがわかる。\n",
    "1. 同じ`np.random`にある`uniform`は`[0,1)`の一様分布から乱数を発生させる関数である。`scipy.stats`の`uniform.rvs`と同じ役割を果たすが，引数の設定が異なる。`low`は最小値，`high`は最大値を指定する。\n",
    "1. `Numpy`の関数である`ones()`は全ての要素が`1`になる$(n\\times 1)$の`array`を生成する。\n",
    "\n",
    "---\n",
    "式（\\*）を使い推定値を計算しよう。まず`Numpy`の関数`stack`を使い`X`を作成する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "X = np.stack([c,x1,x2],axis=1)\n",
    "X.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "ここでの引数`axis=1`は`c`，`x1`，`x2`を縦ベクトルとして横方向につなぎあわせることを指定している。`X`は`Numpy`の`array`であり$(30\\times 3)$の行列となっている。式（\\*）を次のコードで計算する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "bhat = np.linalg.inv((X.T)@X)@(X.T)@y\n",
    "bhat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "左から`b0`，`b1`，`b2`の推定値となる。上のコードを簡単に説明する。\n",
    "* `linalg`は`Numpy`の線形代数のパッケージであり，`inv`は逆行列を計算するための関数\n",
    "* `T`は転置行列の属性\n",
    "* `@`は行列の積\n",
    "\n",
    "`bhat`を計算した１行コードは、定数項以外に2つの説明変数を想定している。実は、定数項以外の説明変数の数が1でも20でも同じコードが使える。ただ、必要なのは`X`を説明変の数に合わせて作成することである。\n",
    "\n",
    "\n",
    "`statsmodels`を使い推定値が同じであることを確認する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "df_check = pd.DataFrame({'Y':y,'X1':x1,'X2':x2})\n",
    "res_check = smf.ols('Y ~ X1 + X2', data=df_check).fit()\n",
    "res_check.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 推定値の標準誤差"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "次に推定値の標準誤差を考えよう。推定値の分散共分散行列は次式で表される。\n",
    "\n",
    "$$\n",
    "\\text{Var}\\left(\\hat{B}\\right)=\\hat{\\sigma}_e^2\\left(X^TX\\right)^{-1}\n",
    "$$\n",
    "\n",
    "$\\sigma_e^2$は次の２つの式から求められる。\n",
    "\n",
    "$$\n",
    "e = Y-X\\hat{B},\\qquad\n",
    "\\text{SSR}=ee^T,\\qquad\n",
    "\\hat{\\sigma}_e^2=\\dfrac{\\text{SSR}}{n-1-k}\n",
    "$$\n",
    "\n",
    "ここで$e$は残差，SSR（Residual Sum of Squares）残差変動の平方和を指しており，$k$は定数項以外の説明変数の数である。係数の推定値の標準誤差は$\\text{Var}\\left(\\hat{B}\\right)$の対角成分の平方根で与えられる。実際に計算してみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "e = y.T-X@bhat.T\n",
    "SSR = e@e.T\n",
    "SSR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "この値は`res_check`の属性`.ssr`と同じことが確認できる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "res_check.ssr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "次に$\\hat{\\sigma}_e^2$を計算しよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "sigma_e2 = SSR/(n-1-2)\n",
    "sigma_e2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "この値も`res_check`の属性`.scale`から取得できる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "res_check.scale"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "次に$\\text{Var}\\left(\\hat{B}\\right)$を計算し変数`bvar`に割り当てる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "bvar = sigma_e2 * np.linalg.inv(X.T@X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`bvar`は3X3の行列になっており，その対角成分はメソッド`.diagonal()`で抽出できる。その平方根が推定値の標準誤差となる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "np.sqrt(bvar.diagonal())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`statsmodels`の結果と比べてみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "res_check.bse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "同じ結果である。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## シミュレーション：不偏性"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "不偏性をシミュレーションで考える。\n",
    "* 標本をランダム抽出するたびに$\\hat{\\beta}_j$, $j=0,1$は異なる値を取るが，標本抽出を繰り返し$\\hat{\\beta}_j$の平均を計算すると真の値と「等しい」ことを示す。\n",
    "\n",
    "シミュレーションでは上で使った同じ回帰式を使う。\n",
    "\n",
    "$$y=\\beta_0+\\beta_1x_1+\\beta_2x_2+u$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "シミュレーションの関数を作成するが，パラメータの推定値の計算には`statsmodels`を使わずに，`Numpy`の関数を使い「手計算」とし高速化パッケージ`Numba`を使う。単回帰分析のシミュレーションと同じように、使い方は簡単でデコレーターと呼ばれる`@njit`（又は`@jit`）を関数の上に加えるだけである。\n",
    "\n",
    "---\n",
    "次の関数の引数：\n",
    "* `n`：標本の大きさ\n",
    "* `N`：シミュレーションの回数（標本数と考えても良い）\n",
    "* `b0`：定数項（デフォルトは`1.0`と設定）\n",
    "* `b1`：スロープ係数（デフォルトは`2.0`と設定）\n",
    "* `b2`：スロープ係数（デフォルトは`3.0`と設定）\n",
    "\n",
    "次の関数の返り値：\n",
    "* `b0`、`b1`、`b2`の`N`個の推定値がそれぞれ格納されている`Numpy`の`array`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "@njit  # 計算の高速化\n",
    "def sim_unbias(n, N, b0=1.0, b1=2.0, b2=3.0):\n",
    "\n",
    "    # N個の0からなるarrayであり、0を推定値と置き換えて格納する\n",
    "    b0hat_arr = np.zeros(N)\n",
    "    b1hat_arr = np.zeros(N)\n",
    "    b2hat_arr = np.zeros(N)\n",
    "    \n",
    "    c = np.ones(n)  # 定数項\n",
    "    \n",
    "    for i in range(N):  # N 回のループ\n",
    "        \n",
    "        # 標準正規分布\n",
    "        u = np.random.normal(loc=0, scale=1.0, size=n)\n",
    "        \n",
    "        # 平均4.0,標準偏差2.0の正規分布\n",
    "        x1 = np.random.normal(loc=4, scale=2.0, size=n)\n",
    "        \n",
    "        # [1,10)の一様分布。`low=`，`high=`を付け加えるとエラーが発生する。\n",
    "        x2= np.random.uniform(1.0, 10.0, size=n)\n",
    "        \n",
    "        y = b0 + b1*x1 + b2*x2 + u  # 被説明変数の生成\n",
    "        \n",
    "        X = np.stack((c,x1,x2),axis=1)  # 説明変数の行列\n",
    "        bhat = np.linalg.inv((X.T)@X)@(X.T)@y  # OLS推定\n",
    "        \n",
    "        b0hat_arr[i] = bhat[0]  # 上のarrayに推定値を代入する\n",
    "        b1hat_arr[i] = bhat[1]\n",
    "        b2hat_arr[i] = bhat[2]\n",
    "        \n",
    "    return b0hat_arr, b1hat_arr, b2hat_arr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "標本の大きさは`30`、シミュレーションの回数は`100_000`に設定し、結果は`b0hat`, `b1hat`, `b2hat`に割り当てる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "b0hat, b1hat, b2hat = sim_unbias(n=30, N=100_000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "推定値の平均を計算してみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "print('b0:', b0hat.mean(),\n",
    "      '\\nb1:', b1hat.mean(),\n",
    "      '\\nb1:', b2hat.mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "平均は真の値と全く同じではないが，十分に近い値である。即ち，不偏性が満たされている。標本数`N`を増加させるとより真の値に近づくことになる。\n",
    "\n",
    "次に、OLS推定量$\\hat{\\beta}_1$の分布を図示しよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "plt.hist(b1hat,bins=30)\n",
    "plt.axvline(x=b1,color='red')\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`b1`の真の値`2.0`を中心に概ね左右対象に分布していることが分かる。\n",
    "\n",
    "次に分布（ヒストグラム）のカーネル密度推定をおこなうために，`scipy.stats`にある`gaussian_kde`を使う。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "x=np.linspace(1.5,2.5,100)  # 図を作成するために1.5から2.5までの横軸の値を設定\n",
    "kde_model=gaussian_kde(b1hat)  # カーネル密度推定を使いb1hatの分布を推定\n",
    "\n",
    "plt.plot(x, kde_model(x))  # 誤差項の分布をプロット\n",
    "plt.axvline(x=b1,color='red')  # 母集団のパラメータ\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "２つの図を重ねる。\n",
    "\n",
    "（注意）ヒストグラムの縦軸は頻度である。一方，カーネル密度推定の場合，曲線の下の面積が１になるように縦軸が設定されている。ヒストグラムの縦軸をカーネル密度推定に合わせるために`density=True`のオプションを加える。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "x=np.linspace(1.5,2.5,100)\n",
    "kde_model=gaussian_kde(b1hat)\n",
    "\n",
    "plt.plot(x, kde_model(x))\n",
    "plt.hist(b1hat,bins=30, density=True)\n",
    "plt.axvline(x=b1,color='red')\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "＜Give it a try!＞\n",
    "* 標本の大きさ`n`を30、100、500の３つのケースでシミュレーションをおこない、それぞれのカーネル密度推定を重ねて１つの図を作成しなさい。\n",
    "* $\\hat{\\beta}_0$と$\\hat{\\beta}_2$の分布も図示しなさい。\n",
    "* `n`を変化させて違いを確かめなさい。\n",
    "* `N`を変化させて違いを確かめなさい。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "tags": [
     "remove-cell"
    ]
   },
   "source": [
    "## シミュレーション：多重共線性"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "(sec:9-vif)=\n",
    "## シミュレーション：多重共線性"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 説明"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "多重共線性が高くなると，パラメータの推定値の精度が低下する（即ち、パラメータの標準誤差$\\text{se}\\left(\\hat{\\beta}_j\\right)$が大きくなる）。この点を確認するためにパラメータの分散の推定値にある\n",
    "\n",
    "$$\n",
    "\\frac{1}{1-R_j^2}\n",
    "$$\n",
    "\n",
    "を利用する。これは\n",
    "**分散拡大因子（variance inflation factor）**\n",
    "と呼ばれる。以下では`vif`と呼ぶ。この値が１０以上になる説明変数どうしを使うと多重共線性の可能性が高いといわれる。\n",
    "\n",
    "例として回帰分析`res_1`を考える。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 方法１：手計算"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "まず説明変数だけから構成される`DataFrame`を作成する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "wage1_vif = wage1.loc[:,['educ','tenure','exper']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`corr()`は相関係数を返すメソッドであり，そこから`.to_numpy()`を使いデータを`Numpy`の`array`として取り出す。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "mc = wage1_vif.corr().to_numpy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "* 上でも使った`linalg.inv()`は、`NumPy`にある`linalg`サブパッケージの関数であり、逆行列を返す。\n",
    "* `diagonal()`は対角成分を返すメソッドであり，それが`vif`である。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "vif_manual = np.linalg.inv(mc).diagonal()\n",
    "vif_manual"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "変数名と一緒に表示するために次のように`Series`として表示する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "pd.Series(vif_manual, index=wage1_vif.columns)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "次に上の計算を関数にまとめる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def my_vif(dataframe):\n",
    "    mc = dataframe.corr().to_numpy()\n",
    "    vif = np.linalg.inv(mc).diagonal()\n",
    "    return pd.Series(vif_manual, index=dataframe.columns)\n",
    "\n",
    "my_vif(wage1_vif)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 方法２：`statsmodels`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`statsmodels`にある`variance_inflation_factor`を使い`vif`を計算することもできる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "必ず**定数項を含む**説明変数だけから構成される`DataFrame`を作成する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "wage1_vif['Intercept'] = 1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "for i in range(len(wage1_vif.columns)-1):  # 定数項は無視するために-1\n",
    "    name = wage1_vif.columns[i]\n",
    "    vif_val = vif(wage1_vif.to_numpy(), i)\n",
    "    print(f'{name : <10}{vif_val}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "定数項は考えなくて良い。\n",
    "\n",
    "（コードの説明）\n",
    "* `f-string`を使っている。\n",
    "* `{name : <10}`は、`name`の値を代入し、`10`という幅の中で左寄せしている。\n",
    "* `{vif_val}`は、`vif_val`の値を代入している。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### シミュレーション"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "シミュレーションをおこない、多重共線性により`OLS`推定量の標準誤差が上昇し推定量の精度が損なわれることを確認する。\n",
    "\n",
    "シミュレーションでは以下の回帰式を使う。\n",
    "\n",
    "$$y=\\beta_0+\\beta_1x_1+\\beta_2x_2+u$$\n",
    "\n",
    "母集団の係数の値を以下のように設定する\n",
    "* $\\beta_0=1.0$\n",
    "* $\\beta_1=2.0$\n",
    "* $\\beta_2=3.0$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "シミュレーションの関数を作成する。\n",
    "* 引数\n",
    "    * `n`：標本の大きさ\n",
    "    * `N`：標本数（ループの回数）\n",
    "    * `m`：２つの説明変数の共分散\n",
    "    * `b0`：定数項（デフォルトは`1.0`と設定）\n",
    "    * `b1`：スロープ係数（デフォルトは`2.0`と設定）\n",
    "    * `b2`：スロープ係数（デフォルトは`3.0`と設定）\n",
    "* 返り値\n",
    "    * $\\hat{\\beta}_i$，$i=0,1,2$の推定値のリスト\n",
    "    \n",
    "（コメント）\n",
    "\n",
    "* `@njit`を使いたいところだが，`Numpy`の`random.multivariate_normal()`が`Numba`に対応していないため`scipy.stats`の `multivariate_normal.rvs()`を使う。`np.random.normal()`を使って二変量正規分布からの値とする方法もあるが，ここでは簡単化を重視する。\n",
    "* `ols`は係数の推定値だけではなく他の多くの統計値も自動的に計算するため一回の計算に比較的に長い時間を要する。計算の速度を少しでも早めるために下の関数の中では`ols`は使わず`Numpy`の関数を使いOLS推定値を計算する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def sim_multi(n, N, m, b0=1.0, b1=2.0, b2=3.0):  # n=標本の大きさ, N=標本数, m=共分散\n",
    "    \n",
    "    # ２つのx1,x2の共分散を設定\n",
    "    rv_mean = [4, 1]  # x1, x2の平均\n",
    "    # x1, x2の共分散行列\n",
    "    rv_cov = [[1.0, m],    # 全ての変数の分散は１（対角成分）\n",
    "              [m, 1.0]]  # Cov(x1,x2)=m\n",
    "    \n",
    "    # 推定値を入れる空のarray\n",
    "    b0hat_arr = np.zeros(N)\n",
    "    b1hat_arr = np.zeros(N)\n",
    "    b2hat_arr = np.zeros(N)\n",
    "    \n",
    "    c = np.ones(n)  # 定数項\n",
    "    \n",
    "    for i in range(N):  # N 回のループ\n",
    "        \n",
    "        # x1, x2の値の抽出\n",
    "        rv = multivariate_normal.rvs(rv_mean, rv_cov, size=n)  # x1, x2,をnセット抽出\n",
    "        x1 = rv[:,0]  # 説明変数\n",
    "        x2 = rv[:,1]  # 説明変数\n",
    "\n",
    "        u = np.random.randn(n)  # 標準正規分布\n",
    "        y = b0 + b1*x1 + b2*x2 + u  # 説明変数\n",
    "        \n",
    "        X = np.stack((c,x1,x2),axis=1)  # 説明変数の行列\n",
    "        bhat = np.linalg.inv((X.T)@X)@(X.T)@y  # OLS推定\n",
    "        \n",
    "        b0hat_arr[i] = bhat[0]  # b0hat_arrへの追加\n",
    "        b1hat_arr[i] = bhat[1]  # b1hat_arrへの追加\n",
    "        b2hat_arr[i] = bhat[2]  # b2hat_arrへの追加\n",
    "\n",
    "    return b0hat_arr, b1hat_arr, b2hat_arr # 返り値の設定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "シミュレーションの開始"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# 多重共線性が弱いケース \n",
    "b0hat_weak, b1hat_weak, b2hat_weak = sim_multi(30, 10_000, m=0.1)\n",
    "\n",
    "# 多重共線性が強いケース \n",
    "b0hat_strong, b1hat_strong, b2hat_strong = sim_multi(30, 10_000, m=0.9)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "$\\hat{\\beta}_1$の分布の図示"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "xx=np.linspace(0.5,3.5,num=100)  # 図を作成するために横軸の値を設定\n",
    "\n",
    "# 多重共線性が弱いケース\n",
    "kde_model_weak=gaussian_kde(b1hat_weak)  # OLS推定量のカーネル密度関数を計算\n",
    "\n",
    "# 多重共線性が強いケース\n",
    "kde_model_strong=gaussian_kde(b1hat_strong)\n",
    "\n",
    "plt.plot(xx, kde_model_weak(xx), 'g-', label='Low VIF')  # OLS推定量の分布プロット\n",
    "plt.plot(xx, kde_model_strong(xx),'r-', label='High VIF')  # 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": [
    "多重共線性が強いと推定値の分布は，真の値（$\\beta_1=$\n",
    "`2.0`）の周辺で低くなり左右に広がっている。推定値の精度が低下することを示している。$\\hat{\\beta}_1$の分散を計算してみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "np.var(b1hat_weak), np.var(b1hat_strong)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "＜Give it a try!＞\n",
    "\n",
    "以下のシミュレーションをおこない，違いは何か確認しなさい。\n",
    "1. $\\hat{\\beta}_0$と$\\hat{\\beta}_2$の分布も図示しなさい。\n",
    "1. `n`を変化させて違いを確かめなさい。\n",
    "1. `N`を変化させて違いを確かめなさい。\n",
    "1. `sim_multi`を書き換えて$\\hat{\\beta}_1$と$\\hat{\\beta}_2$の標準誤差をプロットし問4のそれぞれの平均を計算しなさい。（ヒント：推定結果の属性`bse`を使う）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "問4の答え：$+$を押すと答えが表示される。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true,
    "tags": [
     "hide-cell"
    ]
   },
   "outputs": [],
   "source": [
    "def sim_se(n, N, m, b0=1.0, b1=2.0, b2=3.0):  # n=標本の大きさ, N=標本数, m=共分散\n",
    "    \n",
    "    rv_mean = [4, 1]  # x1, x2の平均\n",
    "    # x1, x2の共分散行列\n",
    "    rv_cov = [[1.0, m],    # 全ての変数の分散は１（対角成分）\n",
    "              [m, 1.0]]  # Cov(x1,x2)=m\n",
    "    \n",
    "    # 推定値を入れるarray\n",
    "    se_x1_arr = np.zeros(N)\n",
    "    se_x2_arr = np.zeros(N)\n",
    "\n",
    "    for i in range(N):  # N 回のループ\n",
    "        \n",
    "        rv = multivariate_normal.rvs(rv_mean, rv_cov, size=n)  # x1, x2,をnセット抽出\n",
    "        x1 = rv[:,0]  # 説明変数\n",
    "        x2 = rv[:,1]  # 説明変数\n",
    "        X = np.stack([np.ones(n),x1,x2],axis=1)  # 説明変数の行列    \n",
    "        \n",
    "        u = np.random.normal(loc=0, scale=1, size=n)  # 誤差項（標準正規分布）\n",
    "        y = b0 + b1*x1 + b2*x2 + u  # 説明変数\n",
    "\n",
    "        bhat = np.linalg.inv((X.T)@X)@(X.T)@y # 係数の推定値\n",
    "        yhat = X@bhat  # yの予測値\n",
    "        uhat = y-yhat  # 残差\n",
    "        rss = (uhat.T)@uhat  # 残差平方和\n",
    "        sigma2 = rss/(n-3)  # 回帰の残差（不偏）分散 \n",
    "        XTX_1 = np.linalg.inv(X.T@X)  # moment matrix\n",
    "        variance = sigma2*np.diag(XTX_1)  # bの分散\n",
    "        se = np.sqrt(variance)  # bの標準誤差\n",
    "        \n",
    "        se_x1_arr[i] = se[1]  # arrayに追加\n",
    "        se_x2_arr[i] = se[2]  # array追加\n",
    "        \n",
    "    return se_x1_arr, se_x2_arr # 返り値の設定\n",
    "\n",
    "\n",
    "# 多重共線性が弱いケース \n",
    "se_x1_weak, se_x1_weak = sim_se(30,1000,m=0.1)\n",
    "\n",
    "# 多重共線性が強いケース \n",
    "se_x1_strong, se_x1_strong = sim_se(30,1000,m=0.9)\n",
    "\n",
    "\n",
    "plt.plot(se_x1_weak, label='Low VIF')\n",
    "plt.plot(se_x1_strong, label='High VIF')\n",
    "plt.title('Standard Error of x1')\n",
    "plt.xlabel('Simulation rounds')\n",
    "plt.legend()\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "問5の答え：$+$を押すと答えが表示される。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true,
    "tags": [
     "hide-cell"
    ]
   },
   "outputs": [],
   "source": [
    "sum(se_x1_weak)/len(se_x1_weak), sum(se_x1_strong)/len(se_x1_strong)\n",
    "\n",
    "# 次のコードでもOK\n",
    "# np.array(se_x1_weak).mean(), np.array(se_x1_strong).mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 図を使ってチェック"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "多重共線性は説明変数間の相関が高いと発生するが，ここでは相関度を図を使ってチェックする。色々な方法があるが，以下では２つ紹介する。\n",
    "1. `pandas`の関数`scatter_matrix`を使う方法\n",
    "1. `seaborn`というパッケージの関数`pairplot`を使う方法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### `pandas`を使う方法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "まずコードの中で`matplotlib`を明示的に導入せずに`pandas`のみを使い図示する方法を紹介する。こちらの方が簡単と感じるかもしれない。散布図を描いてみる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "wage1.plot.scatter('educ','wage')\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`pandas`の`DataFrame`や`Series`には図を作成するメソッド`plot()`が用意されている（裏では`matplotlib`が動いている）。例えば，`wage1`にある`educ`のヒストグラムであれば次のコードで図示することができる。\n",
    "```\n",
    "wage1['educ'].plot(kind='hist')\n",
    "\n",
    "もしくは\n",
    "\n",
    "wage1['educ'].plot.hist()\n",
    "```\n",
    "また`educ`と`wage`の散布図は次のコードで描くことができる。\n",
    "```\n",
    "wage1.plot('educ','wage',kind='scatter')\n",
    "\n",
    "もしくは\n",
    "\n",
    "wage1.plot.scatter('educ','wage')\n",
    "```\n",
    "\n",
    "また１つの`DataFrame`にある複数列データを使い，複数の図を並べることも可能である。興味がある人は[このリンク](https://pandas.pydata.org/pandas-docs/stable/user_guide/visualization.html)を参照しよう。\n",
    "\n",
    "---\n",
    "変数の相関度をチェックするために`pandas.plotting`の`scatter_matrix`を使う。このモジュールは`DataFrame`を引数とする。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "scatter_matrix(wage1)\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "横軸と縦軸のラベルを確認すること。対角線上はそれぞれの変数のヒストグラム，対角線以外の図は縦横のペアの変数の散布図となっており，相関度をある程度目視で確認できる。\n",
    "\n",
    "---\n",
    "次に主なオプションとして２つを紹介する。\n",
    "1. 図の大きさは`figsize=(9, 6)`で指定する。この例では`9`が横幅，`6`が縦幅である。 \n",
    "1. `diagonal='kde'`を指定すると対角線上のヒストグラムをカーネル密度推定に変更できる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "scatter_matrix(wage1, figsize=(9, 6), diagonal='kde')\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "次に相関係数を簡単に計算する方法を紹介する。`DataFrame`のメソッド`corr()`を使うと変数の相関係数を`DataFrame`として返す。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "mat = wage1.corr()\n",
    "mat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "上の図と相関係数の値を見比べて，概ねどのような関係にあるのかを確かめよう。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### `seaborn`を使う方法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`seaborn`というパッケージを使うとより「見栄えがする」図を描くことができる。`seaborn`については[このリンク](https://seaborn.pydata.org)を参照しよう。`seaborn`は`sns`としてインポートするのが慣例である。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "まず上で`DataFrame`のメソッド`corr()`を使い変数の相関係数を計算したが，`seaborn`の`heatmap()`関数を使うと相関係数を色に変換してより見やすい表示となる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "sns.heatmap(mat, vmin=-1, annot=True, cmap='coolwarm')\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "赤はプラス，青はマイナスを示し，色の濃淡は絶対値に連動している。ここで使った３つのオプション（設定しなくても良い）の説明する。\n",
    "* `vmin`：右の縦長の棒は表示される範囲を示すが，その最低値を`-1`にする（デフォルトは自動で設定）\n",
    "    * 最高値を設定する`vmax`もあるが，`.corr()`の対角成分は`1`なので設定する必要はない。\n",
    "* `annot`は相関係数を表示する（デフォルトは`False`）\n",
    "* `cmap`は色を設定する（デフォルトは`None`）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`seaborn`には`matplotlib`の相関度をチェックする`scatter_matrix`関数に対応する`pairplot`があり，より使い勝手が良いと感じるかもしれない。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "sns.pairplot(wage1, height=1.5, aspect=1.3)\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "図の解釈は`pandas`の場合と同じである。\n",
    "\n",
    "`pairplot()`の主な引数：\n",
    "1. それぞれの図（ファセット）のサイズ\n",
    "    * `height`は高さを指定\n",
    "    * 横幅は`aspect`$\\times$`height`で設定する。\n",
    "1. 対角線上のヒストグラムをカーネル密度推定に変更する場合は`diag_kind=kde`と指定する。\n",
    "1. `kind='reg'`を追加すると，散布図に回帰直線が追加され相関度の確認がより簡単になる。\n",
    "\n",
    "`seaborn`の回帰分析に関係するプロットについては[このリンク](https://seaborn.pydata.org/tutorial/regression.html#regression-tutorial)が役立つだろう。"
   ]
  }
 ],
 "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.15.0"
  },
  "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
}