18_Zero_Conditional_Mean.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(chap:18-GM4)=\n",
    "# GM仮定４が満たされない場合"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div name=\"html-admonition\" style=\"font-size: 0.8em\">\n",
    "<input type=\"button\" onclick=\"location.href='https://translate.google.com/translate?hl=&sl=ja&tl=en&u='+window.location;\" value=\"Google translation\" style=\"color:#ffffff;background-color:#008080; height:25px\" onmouseover=\"this.style.background='#99ccff'\" onmouseout=\"this.style.background='#008080'\"/> in English or the language of your choice.\n",
    "</div><br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.formula.api as smf\n",
    "\n",
    "from numba import njit\n",
    "from scipy.stats import multivariate_normal, norm, uniform\n",
    "\n",
    "# 警告メッセージを非表示\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## 説明"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "仮定４：Zero conditional mean $\\text{E}\\left(u|X\\right)=0$\n",
    "\n",
    "* この仮定が満たされない場合，`OLS`推定量は不偏性・一致性が満たされない。\n",
    "* 経済学の実証分析では，この仮定が満たされない場合が多々ある。その主な理由に次の３つが挙げられる。\n",
    "    * 欠落変数\n",
    "    * 測定誤差\n",
    "    * 同時方程式\n",
    "\n",
    "以下ではこの３つについて説明し，シミュレーションで発生するバイアスを確認する。\n",
    "\n",
    "その前に，仮定４が満たされる場合と満たされない場合の違いをシミュレーションを通して直感的に理解する。"
   ]
  },
  {
   "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": [
    "パラメータの真の値を設定する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "b0=1.0\n",
    "b1=1.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "シミュレーションの関数を作成"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def sim_a4(n, m):  # n=標本の大きさ, m=共分散\n",
    "    \n",
    "    rv_mean = [4, 0]  # x, uの平均\n",
    "    rv_cov = [[1.0, m],   # x, uの共分散行列\n",
    "               [m, 0.1]]  # Cov(x,u)=m\n",
    "    \n",
    "    rv = multivariate_normal.rvs(rv_mean, rv_cov, size=n)  # x1, uをnセット抽出\n",
    "    x = rv[:,0]  # 説明変数\n",
    "    u = rv[:,1]  # 誤差項\n",
    "    y = b0 + b1*x + u  # 説明変数\n",
    "\n",
    "    df = pd.DataFrame({'Y':y, 'X':x})  # DataFrameの作成\n",
    "    res = smf.ols(formula='Y ~ X', data=df).fit()  # OLS推定\n",
    "\n",
    "    return x, y, u, res.params.iloc[0], res.params.iloc[1] # 返り値の設定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "ここで重要な役割を果たすのが説明変数と誤差項の共分散を表す`m`である。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 仮定４が満たされる場合：$m=0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "x_0, y_0, u_0, b0hat_0, b1hat_0 = sim_a4(100, m=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "母集団の誤差項$u$と説明変数$x$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "plt.scatter(x_0,u_0)\n",
    "plt.axhline(0)\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "共分散を計算してみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "np.corrcoef(x_0,u_0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "共分散は非常に小さな数字であり，特定な傾向がない。次に，母集団回帰直線とOLS回帰直線を比べてみる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "xx = np.linspace(min(x_0),max(x_0),100)\n",
    "\n",
    "plt.scatter(x_0,y_0)\n",
    "plt.plot(xx,b0hat_0+b1hat_0*xx, 'k', label='OLS')\n",
    "plt.plot(xx,b0+b1*xx,'r', label='Pop Regression Line')\n",
    "plt.legend()\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "推定量が不偏性を満たしているため両方は殆ど同じである。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 仮定４が満たされない場合：$m\\neq 0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "x_1, y_1, u_1, b0hat_1, b1hat_1 = sim_a4(100, m=0.25)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "母集団の誤差項$u$と説明変数$x$を図示しよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "plt.scatter(x_1,u_1)\n",
    "plt.axhline(0)\n",
    "pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "np.corrcoef(x_1,u_1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "強い正の共分散が確認できる。母集団回帰線とOLS回帰直線を比べる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "xx = np.linspace(min(x_1),max(x_1),100)\n",
    "\n",
    "plt.scatter(x_1,y_1)\n",
    "plt.plot(xx,b0hat_1+b1hat_1*xx, 'k', label='OLS')\n",
    "plt.plot(xx,b0+b1*xx, 'r', label='Pop Regression Line')\n",
    "plt.legend()\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "明らかに異なる。GM仮定４が満たされず，推定量の不偏性が満たされないためである。この場合，一致性も満たされない。標本の大きさ`n`を`100000`にして確かめてみる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "x_1, y_1, u_1, b0hat_1, b1hat_1 = sim_a4(100_000, m=0.25)\n",
    "\n",
    "xx = np.linspace(min(x_1),max(x_1),100)\n",
    "\n",
    "plt.scatter(x_1,y_1, alpha=0.1)  # 下の説明を参照\n",
    "plt.plot(xx,b0hat_1+b1hat_1*xx, 'k', label='OLS')\n",
    "plt.plot(xx,b0+b1*xx,'r', label='Pop Regression Line')\n",
    "plt.legend()\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "上から３行目に`alpha=0.1`とあるが，散布図のマーカーの透明度をしている引数であり`0~1`の値をとる。デフォルトは`1`であり，`0`は透明になる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## 欠落変数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 欠落変数バイアス"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "母集団のモデルは以下で与えられるとしよう。\n",
    "\n",
    "$$y=\\beta_0+\\beta_1 x_1+\\beta_2x_2+\\beta_3x_3+u\\qquad\\quad(\\text{母集団モデル})$$\n",
    "\n",
    "実際に推定されるモデル（$x_3$が欠落）：\n",
    "\n",
    "$$y=\\gamma_0+\\gamma_1 x_1+\\gamma_2x_2+\\eta\\qquad\\quad(\\text{推定モデル})$$\n",
    "\n",
    "更に，$x_1$と$x_3$には以下の関係があると仮定する。\n",
    "\n",
    "$$x_3=\\delta_0+\\delta_1 x_1+\\epsilon\\qquad\\quad(\\text{仮定})$$\n",
    "\n",
    "この式を母集団モデルに代入すると：\n",
    "\n",
    "$$y=(\\beta_0+\\beta_3\\delta_0)+(\\beta_1+\\beta_3\\delta_1)x_1+\\beta_2x_2+\\tilde{\\eta}$$\n",
    "\n",
    "ここで，$\\tilde{\\eta}\\equiv\\beta_3\\epsilon+u$。推定モデルはこの関係に基づいて係数を計算することになる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "**（結果１）**\n",
    "\n",
    "推定モデルを使うと\n",
    "* $\\hat{\\gamma}_0$は$\\beta_0+\\beta_3\\delta_0$の推定値\n",
    "* $\\hat{\\gamma}_1$は$\\beta_1+\\beta_3\\delta_1$の推定値\n",
    "* $\\hat{\\gamma}_2$は$\\beta_2$の推定値\n",
    "\n",
    "となり$\\gamma_0$と$\\gamma_1$の推定値にバイアスが発生する。一方，$\\hat{\\gamma}_2$にはバイアスは発生しない。欠落変数と無相関の変数（定数以外）にはバイアスは発生しないためである。\n",
    "\n",
    "**（結果２）**\n",
    "* 欠落変数が回帰式の説明変数と無相関であれば，バイアスは発生しない！"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "**（理由）**\n",
    "\n",
    "GM仮定４が満たされていないため。母集団モデル，推定モデルと仮定の式から\n",
    "\n",
    "$$\n",
    "\\eta=\\beta_3x_3+u=\\beta_3(\\delta_0+\\delta_1x_1+\\epsilon)+u\n",
    "\\quad\\Rightarrow\\quad\n",
    "\\tilde{\\eta}=\\eta-\\beta_3(\\delta_0+\\delta_1x_1)\n",
    "$$\n",
    "\n",
    "これから$x_1$と$\\tilde{\\eta}$は相関することがわかる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### シミュレーション１：推定値の計算"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "One-shotのシミュレーションをおこないバイアスの発生を確認する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# 母集団のパラメータ\n",
    "b0 = 1\n",
    "b1 = 1\n",
    "b2 = 1\n",
    "b3 = 1\n",
    "\n",
    "# 標本数\n",
    "n = 1000\n",
    "\n",
    "# 説明変数　x1, x2\n",
    "x1 = np.random.normal(4, 1, size=n)\n",
    "x2 = np.random.uniform(1, 10, size=n)\n",
    "\n",
    "# 説明変数　x3\n",
    "e = np.random.normal(size=n)\n",
    "m = 2.0  # x1との相関性を捉える\n",
    "x3 = 1 + m*x1 + e\n",
    "\n",
    "# 母集団の誤差項\n",
    "u = np.random.normal(size=n)\n",
    "y = b0 + b1*x1 + b2*x2 + b3*x3 + u\n",
    "\n",
    "# 標本のDataFrame\n",
    "df_sample = pd.DataFrame({'Y':y, 'X1':x1, 'X2':x2, 'X3':x3})\n",
    "\n",
    "# 推定\n",
    "formula_omitted = 'Y ~ X1 + X2'\n",
    "res_omitted = smf.ols(formula_omitted, data=df_sample).fit()\n",
    "res_omitted.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`X1`のパラメータの値は真の値と明らかに異なる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### シミュレーション２：推定値の分布"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`OLS`推定量の分布を確認する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# 母集団のパラメータの値\n",
    "b0 = 1\n",
    "b1 = 1\n",
    "b2 = 1\n",
    "b3 = 1\n",
    "\n",
    "# シミュレーションの回数\n",
    "N = 10_000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "シミュレーションの関数を設定する。\n",
    "* 引数：\n",
    "    * `n`：標本の大きさ\n",
    "    * `m`：`x1`と`x3`の相関を捉えるパラメータ\n",
    "* 戻り値：推定値のリスト"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "@njit\n",
    "def sim_omitted(n,m):\n",
    "\n",
    "    # 空のリスト\n",
    "    b0hat_arr = np.zeros(N)\n",
    "    b1hat_arr = np.zeros(N)\n",
    "    b2hat_arr = np.zeros(N)\n",
    "\n",
    "    # 説明変数\n",
    "    x1 = np.random.normal(loc=4, scale=1, size=n)\n",
    "    x2 = np.random.uniform(1, 10, size=n)\n",
    "    e = np.random.normal(loc=0, scale=1, size=n)\n",
    "    x3 = 1 + m*x1 + e\n",
    "    \n",
    "    c = np.ones(n)  # 定数項\n",
    "    \n",
    "    for i in range(N):  # N回のループ\n",
    "        # 母集団\n",
    "        u = np.random.normal(loc=0, scale=1, size=n)\n",
    "        y = b0 + b1*x1 + b2*x2 + b3*x3 + 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]\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": [
    "シミュレーションの開始"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "b0hat, b1hat, b2hat = sim_omitted(1000,m=2.0)"
   ]
  },
  {
   "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=30)\n",
    "plt.axvline(x=b0, color='red')\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=30)\n",
    "plt.axvline(x=b1, color='red')\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "$\\hat{\\beta}_2$の分布"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "plt.hist(b2hat, bins=30)\n",
    "plt.axvline(x=b2, color='red')\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## 測定誤差"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 測定誤差によるバイアス"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "次の母集団回帰式を考えよう。\n",
    "\n",
    "$$y=\\beta_0+\\beta_1 x^* + \\eta$$\n",
    "\n",
    "ここで，被説明変数$y$は正確に測定できるが，説明変数$x^*$は以下の式に従って測定される仮定する。\n",
    "\n",
    "$$x=x^*+e$$\n",
    "\n",
    "* $x$：測定値\n",
    "* $e$：測定誤差\n",
    "\n",
    "（仮定）\n",
    "\n",
    "* 測定誤差$e$は真の値と無関係。即ち，$\\text{Cov}(x^*,e)=0$\n",
    "\n",
    "（結果）\n",
    "\n",
    "* 次式をOLS推定する場合，$\\hat{\\beta}_1$は不偏性・一致性を満たさない。\n",
    "\n",
    "$$y=\\beta_0+\\beta_1 x + u,\\qquad u=\\eta-\\beta e$$\n",
    "\n",
    "（理由）\n",
    "\n",
    "仮定４：$\\text{Cov}(x,u)=0$が満たされない。\n",
    "\n",
    "$$\n",
    "\\text{Cov}(x,u)=\\text{E}[xu]\n",
    "=\\text{E}\\left[(x^*+e)(\\eta-\\beta e)\\right]\n",
    "=-\\beta\\cdot\\text{E}(e^2)>0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### シミュレーション１：推定値の計算"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "One-shotのシミュレーションをおこないバイアスの発生を確認する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# 標本の大きさ\n",
    "n = 100\n",
    "\n",
    "# 母集団回帰式\n",
    "b0 = 1.0\n",
    "b1 = 1.0\n",
    "x_pop = np.random.uniform(1,10,size=n)  # 母集団の説明変数\n",
    "u = np.random.normal(scale=1, size=n)  # 母集団の誤差項\n",
    "y = b0 + b1*x_pop + u  # 母集団回帰式\n",
    "\n",
    "# 測定誤差\n",
    "error_sd = 3  # 測定誤差の標準偏差\n",
    "x = x_pop + np.random.normal(scale=error_sd, size=n)  # 測定誤差\n",
    "\n",
    "# OLS推定\n",
    "df_sim = pd.DataFrame({'Y':y, 'X':x})  # DataFrameの作成\n",
    "res_ols = smf.ols('Y ~ X', data=df_sim).fit()  # OLS推定\n",
    "\n",
    "res_ols.params  # OLS推定量"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### シミュレーション２：推定値の分布"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`OLS`推定量の分布を確認する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# 真のパラメータ\n",
    "b0 = 1.0\n",
    "b1 = 1.0\n",
    "\n",
    "# シミュレーションの回数\n",
    "N = 100_000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "@njit\n",
    "def sim_measure(n):\n",
    "    \n",
    "    b0hat_arr = np.zeros(N)\n",
    "    b1hat_arr = np.zeros(N)\n",
    "\n",
    "    x_pop = np.random.uniform(1,10,size=n)  # 母集団の説明変数\n",
    "    \n",
    "    c = np.ones(n)  # 定数項\n",
    "        \n",
    "    # 説明変数の測定誤差\n",
    "    error_sd = 3  # 測定誤差の標準偏差\n",
    "    x = x_pop + np.random.normal(loc=0, scale=error_sd, size=n)  # 測定誤差\n",
    "\n",
    "    for i in range(N):\n",
    "        u = np.random.normal(loc=0, scale=1, size=n)  # 母集団の誤差項\n",
    "        y = b0 + b1*x_pop + u  # 母集団回帰式\n",
    "        \n",
    "        X = np.stack((c,x), axis=1)  # 説明変数の行列\n",
    "        bhat = np.linalg.inv((X.T)@X)@(X.T)@y  # OLS推定\n",
    "        \n",
    "        b0hat_arr[i] = bhat[0]\n",
    "        b1hat_arr[i] = bhat[1]\n",
    "\n",
    "    return b0hat_arr, b1hat_arr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "シミュレーションの開始"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "b0hat, b1hat = sim_measure(100)"
   ]
  },
  {
   "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=30)\n",
    "plt.axvline(x=b0,color='red')\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=30)\n",
    "plt.axvline(x=b1,color='red')\n",
    "pass"
   ]
  },
  {
   "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": [
    "同時方程式モデルとは，均衡メカニズムなどを通して複数の内生変数が複数の式によって同時決定されるモデルである。例として労働の需給モデルを考えよう。均衡では需要量（$L_d$）と供給量（$L_s$）は等しくなり（$L=L_d=L_s$），需要と供給はそれぞれ均衡賃金（$W$）に依存する。\n",
    "* 労働供給関数\n",
    "\n",
    "    $$ L = s_0+s_1 W + s_2 X_s + u_s\\qquad\\qquad\\qquad\\text{(式１)}$$\n",
    "\n",
    "    * $s_1>0$\n",
    "    * $X_s=$ 供給の「その他」の決定要因（例えば，所得水準）\n",
    "    * $u_s=$ 供給の誤差項\n",
    "\n",
    "* 労働需要関数\n",
    "\n",
    "    $$ W = d_0+d_1 L + d_2 X_d + u_d\\qquad\\qquad\\qquad\\text{(式２)}$$\n",
    "        \n",
    "    * $d_1<0$\n",
    "    * $X_d=$ 需要の「その他」の決定要因（例えば，教育水準）\n",
    "    * $u_d=$ 需要の誤差項\n",
    "\n",
    "（相関性の仮定）\n",
    "* $\\text{Cov}(X_s,u_s)=\\text{Cov}(X_d,u_d)=\\text{Cov}(u_s,u_d)=0$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "上の式は**構造方程式（structural equations）**\n",
    "と呼ばれる。これらの式を直接推定するとOLS推定量にはバイアスが発生する。理由は，供給関数では$\\text{Cov}(W,u_s)\\neq 0$となり需要関数でも$\\text{Cov}(L,u_d)\\neq 0$となるためである。即ち，仮定４が満たされないのである。この結果は簡単に示すことができる。２つの式を同時方程式として$L$と$W$の解を計算すると次式を得る。\n",
    "\n",
    "$$L=\\alpha_0+\\alpha_1X_d+\\alpha_2X_s+\\frac{s_1u_d+u_s}{1-d_1s_1}\\qquad\\qquad\\qquad\\text{(式３)}$$\n",
    "\n",
    "$$W=\\beta_0+\\beta_1X_d+\\beta_2X_s+\\frac{d_1u_s+u_d}{1-d_1s_1}\\qquad\\qquad\\qquad\\text{(式４)}$$\n",
    "\n",
    "ここで$\\alpha_i$,$\\beta_i$,$i=0,1,2$は$d_i$,$s_i$,$i=0,1,2$の非線形な形をとることになる。\n",
    "このような$L$と$W$の解の式は**誘導型方程式（reduced-form equations）**\n",
    "と呼ばれるが，この式から次の点は明らかである。\n",
    "* （式３）：$L$は$u_d$に依存し，$L$と$u_d$は相関する　$\\Rightarrow$　$\\text{Cov}(L,u_d)\\neq 0$\n",
    "* （式４）：$W$も$u_s$に依存し，$W$と$u_s$は相関する　$\\Rightarrow$　$\\text{Cov}(W,u_s)\\neq 0$\n",
    "\n",
    "では誘導型方程式を推定すればどうなるのか？相関性の仮定のもと$\\hat{\\alpha}_i$,$\\hat{\\beta}_i$,$i=0,1,2$は不偏性・一致性を保持する。しかし問題は，誘導型方程式を満たす構造方程式は無数にあり，その推定値から$d_i$,$s_i$,$i=0,1,2$を復元できないのである。従って，それらのパラメータの値を推定するためには(式１)と(式２)を推定する必要がある。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### シミュレーション１：推定値の計算"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "（目的）\n",
    "\n",
    "One-shotのシミュレーションをおこない構造方程式を直接推定すると推定値にバイアスが発生することを確認する。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "シミュレーションの簡単化のために(式２)の中で$s_2=0$を仮定する。これにより，(式３)と(式４)のパラメータは以下で与えられる。\n",
    "\n",
    "$$\n",
    "\\alpha_0\\equiv\\frac{s_0+s_1d_0}{1-s_1d_1},\\quad\n",
    "\\alpha_1\\equiv\\frac{s_1d_2}{1-s_1d_1},\\quad\n",
    "\\alpha_2=0\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\beta_0\\equiv\\frac{d_0+d_1s_0}{1-s_1d_1},\\quad\n",
    "\\beta_1\\equiv d_2,\\quad\n",
    "\\beta_2=0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "（シミュレーション・プラン）\n",
    "* $s_0$と$d_i$, $i=0,1,2$の値を設定し，(式３)と(式４)に従う抽出データを生成する（標本の大きさ$=n$）\n",
    "* この標本データを使い(式２)を推定し，推定量のバイアスを確認する。\n",
    "\n",
    "母集団のパラメータを次のように設定する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "s0 = 1.0\n",
    "s1 = 0.5\n",
    "\n",
    "d0 = 1.5\n",
    "d1 = -1\n",
    "d2 = 2.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "標本の大きさ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "n = 10_000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "抽出された説明変数と被説明変数（誤差項）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "xd = np.random.uniform(1,10,size=n)  # 説明変数\n",
    "ud = np.random.normal(loc=0, scale=1, size=n)  # 母集団の誤差項\n",
    "us = np.random.normal(loc=0, scale=2, size=n)  # 母集団の誤差項"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "母集団の労働供給曲線"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "l = (s0+s1*d0)/(1-s1*d1) + (s1*d2/(1-s1*d1))*xd + (s1*ud+us)/(1-s1*d1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "母集団の労働需要曲線"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "w = (d0+d1*s0)/(1-s1*d1) + d2*xd + (d1*us+ud)/(1-s1*d1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "労働需要曲線を使い賃金をOLS推定"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "df_simul = pd.DataFrame({'W':w, 'L':l, 'Xd':xd})  # DataFrameの作成\n",
    "\n",
    "res_ols = smf.ols('W ~ L + Xd', data=df_simul).fit()  # OLS推定\n",
    "\n",
    "res_ols.params  # OLS推定量"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "分かりやすく表示する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "print(f'd0の真の値：{d0}\\t\\td0の推定量：{res_ols.params.iloc[0]}')\n",
    "print(f'd1の真の値：{d1}\\t\\td1の推定量：{res_ols.params.iloc[1]}')\n",
    "print(f'd2の真の値：{d2}\\t\\td2の推定量：{res_ols.params.iloc[2]}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "標本の大きさが非常に大きくてもバイアスは消えない。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### シミュレーション２：推定値の分布"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`OLS`推定量の分布を確認する。シミュレーション１と同じパラメータの真の値を使う。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "s0 = 1.0\n",
    "s1 = 0.5\n",
    "\n",
    "d0 = 1.5\n",
    "d1 = -1\n",
    "d2 = 2.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "シミュレーションの回数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "N = 100_000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "シミュレーションの関数を設定する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "@njit\n",
    "def sim_simul(n):\n",
    "    \n",
    "    b0hat_arr = np.zeros(N)\n",
    "    b1hat_arr = np.zeros(N)\n",
    "    b2hat_arr = np.zeros(N)\n",
    "\n",
    "    xd = np.random.uniform(1, 10, size=n)  # 説明変数\n",
    "    \n",
    "    c = np.ones(n)  # 定数項\n",
    "    \n",
    "    for i in range(N):\n",
    "        # 母集団の誤差項\n",
    "        ud = np.random.normal(loc=0, scale=1, size=n)\n",
    "        \n",
    "        # 母集団の誤差項\n",
    "        us = np.random.normal(loc=0, scale=2, size=n)\n",
    "        \n",
    "        # 母集団の労働供給曲線\n",
    "        l = (s0+s1*d0)/(1-s1*d1) + (s1*d2/(1-s1*d1))*xd + (s1*ud+us)/(1-s1*d1)\n",
    "        \n",
    "        # 母集団の労働需要曲線\n",
    "        w = (d0+d1*s0)/(1-s1*d1) + d2*xd + (d1*us+ud)/(1-s1*d1)\n",
    "\n",
    "        # 説明変数の行列\n",
    "        X = np.stack((c,l,xd), axis=1)\n",
    "        \n",
    "        # 労働需要曲線を使い賃金をOLS推定\n",
    "        bhat = np.linalg.inv((X.T)@X)@(X.T)@w  # OLS推定\n",
    "        \n",
    "        b0hat_arr[i] = bhat[0]\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": [
    "シミュレーションの開始"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "b0hat, b1hat, b2hat = sim_simul(100)"
   ]
  },
  {
   "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=30)\n",
    "plt.axvline(x=d0,color='red')\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=30)\n",
    "plt.axvline(x=d1,color='red')\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "$\\hat{\\beta}_2$の分布"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "plt.hist(b2hat,bins=30)\n",
    "plt.axvline(x=d2,color='red')\n",
    "pass"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.9"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}