14_Hetero.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(chap:14-hetero)=\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 lmdiag\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.formula.api as smf\n",
    "import wooldridge\n",
    "\n",
    "from seaborn import residplot\n",
    "from statsmodels.stats.api import het_breuschpagan, het_white\n",
    "from statsmodels.stats.outliers_influence import reset_ramsey\n",
    "\n",
    "# 警告メッセージを非表示\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 説明"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**＜仮定５（均一分散; Homogeneity）が満たされない場合＞**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "仮定５の均一分散の下では説明変数は誤差項の分散に影響を与えない。即ち，\n",
    "\n",
    "$$\\text{Var}\\left(u|x\\right)=\\sigma^2$$\n",
    "\n",
    "この仮定に関連して次の点について留意する必要がある。\n",
    "* 均一分散が満たされない場合でも，\n",
    "    * 仮定１〜４のもとでOLS推定量 $\\hat{\\beta}_j$の不遍性と一致性は成立する。\n",
    "    * $R^2$の解釈は変わらない。\n",
    "* しかし，均一分散が満たされない場合，\n",
    "    * OLS推定量の標準偏差の推定量である標準誤差$\\text{se}\\left(\\hat{\\beta}_j\\right)$は無効となる。従って，$t$検定と$F$検定も無効になる。\n",
    "    * 大標本特性（漸近的特性）も成立しない。従って，大標本であっても$t$検定と$F$検定も無効になる。\n",
    "\n",
    "仮説を検証するということを目的とすると，検定が無効というのは致命的な問題である。特に，不均一分散（Heteroskedasticity）の問題は，横断面データを使うと頻繁に出てくる問題である。ではどのように対応すればよいのか。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "**不均一分散頑健的推定（Heteroskedasticity-Robust Inference）**\n",
    "\n",
    "この手法を使うと，OLS推定量の標準誤差が調整され**未知の**不均一分散であっても，$t$検定，$F$検\n",
    "定が有効になるというものである。\n",
    "\n",
    "（理由）均一分散であっても不均一分散であっても，$n\\rightarrow\\infty$の場合，不均一分散頑健的推定の$t$($F$)値は$t$($F$)分布に従う。言い換えると，標本の大きさが十分に大きければ，$t$($F$)値の分布は$t$($F$)分布で近似できるということである。\n",
    "\n",
    "更なる利点は，通常のOLS推定の後に標準誤差の調整が施され，計算は`statsmodels`を使うと簡単におこなうことが可能である。\n",
    "\n",
    "（注意）大標本でのみ有効。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "不均一分散頑健的推定では，OLS推定の共分散行列（covariance matrix）と呼ばれる箇所を調整し，OLS推定量の標準誤差を修正する。その調整方法が複数提案されていおり，`statsmodels`では以下の種類に対応している。\n",
    "* `HC0`: White (1980)の不均一分散頑健共分散行列推定\n",
    "* `HC1`: MacKinnon and White (1985)の不均一分散頑健共分散行列推定v1\n",
    "* `HC2`: MacKinnon and White (1985)の不均一分散頑健共分散行列推定v2\n",
    "* `HC3`: MacKinnon and White (1985)の不均一分散頑健共分散行列推定v3\n",
    "\n",
    "ここで`HC`は`H`eteroskedasticity-`C`onsistent Covariance Matrix Estimatorsの`H`と`C`。\n",
    "\n",
    "不均一分散頑健共分散行列推定を使い計算した$t$値や$F$値を**頑健的$t$値**，**頑健的$F$値**と呼ぶことにする。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "````{note}\n",
    "`HC0`などは不均一分散に対応する推定方法である。一方で，時系列分析では残差に不均一分散と自己相関の両方が存在する場合がある。この両方の問題に頑健的な推定を不均一分散・自己相関頑健推定（Heteroskedasticity-Autocorrelation Robust Inference）と呼ぶ。次の引数を使うことにより，不均一分散・自己相関頑健標準誤差を計算することができる。\n",
    "```\n",
    "cov_type='HAC', cov_kwds={'maxlags': n}, use_t=True\n",
    "```\n",
    "* `cov_type='HAC'`：不均一分散・自己相関頑健推定の指定\n",
    "* `cov_kwds={'maxlags': n}`（必須）：`n`は残差のラグを示し，`n`の値に基づきパラメータの推定値の標準誤差を調整する。`n=10`とすれば`10`までのラグを考慮し計算する。標本の大きさより小さな値を設定する必要がある。`for`ループを使い`n`を変化させて、標準誤差がどう変化するかを確認するのも良いだろう。また、`n`を決める一つのルールとして次を紹介しよう。\n",
    "    * Newey and West (1987)が提案したルール：`N`を標本の大きさとすると、次の値の整数部分を`n`に使う。\n",
    "\n",
    "$$\n",
    "4\\left(\\frac{N}{100}\\right)^{2/9}の整数部分\n",
    "$$\n",
    "\n",
    "* `use_t=True`（必須ではない）：`t`検定を使うことを指定する。\n",
    "\n",
    "ちなみに，`HC`は`H`eteroskedasticity `C`onsistentの略であり，`HAC`は`H`eteroskedasticity-`A`utocorrelation `C`onsistentの略である。\n",
    "````"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "OLS推定量の不均一分散頑健標準偏差が簡単に計算できるのであれば，通常の標準偏差を使う必要はないのではないか，という疑問が生じる。この問を関して以下の点を考える必要がある。\n",
    "* 通常の標準偏差を使う利点\n",
    "    * 均一分散の場合（仮定1〜6（CLR仮定）），標本の大きさの大小に関わらず，$t$ ($F$)値の分布は**厳密に** $t$ ($F$)分布に従う。\n",
    "* 不均一分散頑健標準偏差\n",
    "    * 小標本の場合\n",
    "        * 頑健的 $t$ ($F$)値の分布は必ずしも $t$ ($F$)分布に従うわけではない。その場合，$t$ ($F$)検定は無効となる。\n",
    "    * 大標本の場合\n",
    "        * $t$ ($F$)値の分布は $t$ ($F$)分布で**近似**され，$t$ ($F$)検定は有効である。\n",
    "        * この結果は仮定１〜６（CLR仮定）のもとでも同じ。\n",
    "\n",
    "従って，標本の大きさが「大標本」と判断できる場合（例えば，$n\\geq 1000$）以外は通常の標準偏差と不均一分散頑健標準偏差の両方を表示することを勧める。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 頑健的$t$検定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`wooldridge`パッケージのデータセット`gpa3`を使い説明する。この例では大学のGPAと高校の成績や性別，人種などがどのような関係にあるかを探る。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "output_scroll"
    ]
   },
   "outputs": [],
   "source": [
    "gpa3 = wooldridge.data('gpa3').query('spring == 1')  # 春学期だけを抽出\n",
    "wooldridge.data('gpa3', description=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`gpa2`に一部の変数の説明が続いている。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "output_scroll"
    ]
   },
   "outputs": [],
   "source": [
    "wooldridge.data('gpa2', description=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### OLS推定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "被説明変数：\n",
    "* `cumgpa`：累積GPA\n",
    "\n",
    "説明変数\n",
    "* `sat`：SATの成績\n",
    "* `hsperc`：高校の成績の％点（上位から）\n",
    "* `tothrs`：データ抽出時から学期までの時間？（`gpa3`の定義）\n",
    "* `female`：女性ダミー変数（女性＝`1`）\n",
    "* `black`：人種ダミー変数（黒人＝`1`）\n",
    "* `white`：人種ダミー変数（白人＝`1`）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "form_ols = 'cumgpa ~ sat + hsperc + tothrs + female + black + white'\n",
    "\n",
    "mod_ols = smf.ols(form_ols, data=gpa3)\n",
    "res_ols = mod_ols.fit()\n",
    "\n",
    "print(res_ols.summary().tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 不均一分散頑健推定：方法１"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上のOLSの結果を使い頑健$t$値を計算するために，`res_ols`のメソッド`.get_robustcov_results()`を使う。\n",
    "* オプション`cov_type`は頑健性の計算法の指定（デフォルトは`CH1`）。\n",
    "* オプション`use_t`は$t$検定を指定（デフォルトは`None`で「自動」に決められる）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_robust = res_ols.get_robustcov_results(cov_type='HC3', use_t=True)\n",
    "\n",
    "print(res_robust.summary().tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. OLS推定量`coef`の値は同じであり，必ずそうなる。不均一分散頑健推定は，表の中で標準誤差，$t$値，$p$値，信頼区間に影響を与える。\n",
    "2. 標準誤差`std err`を比べると，帰無仮説$\\hat{\\beta}_j=0$の棄却判断を覆すほど大きく変わる変数はない。これは不均一分散がそれほど大きな問題ではないことを示唆している。この点を確かめるために，`res_ols`の誤差項を図示してみる。\n",
    "\n",
    "誤差項を図示する方法として２つを紹介する。一つ目は，`res_ols`の属性`.resid`を使う。`res_ols.resid`は`Pandas`の`Series`（シリーズ）なので，そのメソッド`plot()`を使い図示する。`style`はマーカーを指定するオプション。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_ols.resid.plot(style='o')\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "２つ目の方法は`plt.plot()`を使う。オップション`'o'`はマーカの指定である。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(res_ols.resid, 'o')\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "３つ目は`scatter()`を使う。`.index`は`res_ols.resd`の属性でインデックスを示す。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(res_ols.resid.index, res_ols.resid)\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### 不均一分散頑健推定：方法２"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "OLS推定をする際，`fit()`の関数に`.get_robustcov_results()`で使った同じオプションを追加すると頑健$t$値などを直接出力できる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "res_HC3 = smf.ols(form_ols, data=gpa3).fit(cov_type='HC3', use_t=True)\n",
    "\n",
    "print(res_HC3.summary().tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## 頑健的$F$検定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "同じデータ`gpa3`を使い，黒人ダミーと白人ダミーの係数は両方とも`0`という仮説を検定する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "hypotheses = 'black = white = 0'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "まず通常の$F$検定を考える。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "f_test_ols = res_ols.f_test(hypotheses)\n",
    "\n",
    "f_test_ols.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "返り値（左から）\n",
    "* `F statistic`：$F$統計量\n",
    "* `F p-value`：$F$の$p$値\n",
    "* `df_denom`：分母の自由度\n",
    "* `df_num`：分子の自由度\n",
    "\n",
    "次に頑健$F$検定の方法を説明する。上の不均一分散頑健推定の方法２で使った`f_test_HC3`を使う。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "f_test_HC3 = res_HC3.f_test(hypotheses)\n",
    "\n",
    "f_test_HC3.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "\n",
    "$t$検定の場合と同じように，大きく変わる結果につながってはない。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## 均一分散の検定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "均一分散の場合 $t$($F$)値は厳密に$t$($F$)分散に従う。それが故に，均一分散が好まれる理由である。ここでは均一分散の検定について考える。帰無仮説と対立仮説は以下となる。\n",
    "\n",
    "$\\text{H}_0$：誤差項は均一分散\n",
    "\n",
    "$\\text{H}_A$：帰無仮説は成立しない\n",
    "\n",
    "２つの検定方法を考える。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### ブルーシュ・ペーガン（Breusch-Pagan）検定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "データ`hprice1`を使って，住宅価格の決定要因を検討する。ここで考える均一分散の検定にBreusch-Pagan検定と呼ばれるもので，`statsmodels`のサブパッケージ`stats`の関数`het_breuschpagan`を使う。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "hprice1 = wooldridge.data('hprice1')\n",
    "wooldridge.data('hprice1', description=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "以下で使う変数について。\n",
    "\n",
    "被説明変数\n",
    "* `price`：住宅価格（単位：1000ドル）\n",
    "\n",
    "説明変数\n",
    "* `lotsize`：土地面積（単位：平方フィート）\n",
    "* `sqrft`：家の面積（単位：平方フィート）\n",
    "* `bdrms`：寝室の数"
   ]
  },
  {
   "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": [
    "form_h = 'price ~ lotsize + sqrft + bdrms'\n",
    "\n",
    "res_h = smf.ols(form_h, data=hprice1).fit()\n",
    "\n",
    "print(res_h.summary().tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "この結果に対してBreusch-Pagan検定をおこなう。`het_breuschpagan()`の引数について：\n",
    "* 第１引き数：OLSの結果`res_h`の属性`.resid`で誤差項の値\n",
    "* 第２引き数：OLSの結果`res_h`の属性`.model`の属性`exog`を使い定数項を含む説明変数の値"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "het_breuschpagan(res_h.resid, res_h.model.exog)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "返り値（上から）\n",
    "* `LM statistic`：$LM$統計量\n",
    "* `LM p-value`：$LM$の$p$値\n",
    "* `F statistics`：$F$統計量\n",
    "* `F p-value`：$F$の$p$値\n",
    "\n",
    "$LM$検定とはLagrange Multiplier検定のことで，大標本の場合に仮定１〜４（GM仮説）のもとで成立する。一般にBreusch-Pagan検定は$LM$統計量を使ったものを指すが，$F$統計量としても計算できる。\n",
    "\n",
    "5%の有意水準で帰無仮説（$\\text{H}_0$：誤差項は均一分散）を棄却でき，不均一分散の可能性が高い。対処法として変数の変換が挙げられる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### 対数変換"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`bdrms`以外を対数変換する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "form_h_log = 'np.log(price) ~ np.log(lotsize) + np.log(sqrft) + bdrms'\n",
    "\n",
    "res_h_log = smf.ols(form_h_log, data=hprice1).fit()\n",
    "\n",
    "print(res_h_log.summary().tables[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "het_breuschpagan(res_h_log.resid, res_h_log.model.exog)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "5%の有意水準でお帰無仮説を棄却できない。即ち，対立仮説（$\\text{H}_A$：帰無仮説は成立しない）を採択し，均一分散の可能性が高い。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### ホワイト（White）検定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "この検定はOLS推定量の標準誤差と検定統計量を無効にする不均一分散を主な対象としており，Breusch-Pagan検定よりもより一般的な式に基づいて検定をおこなう。`statsmodels`のサブパッケージ`stats`の関数`het_white`を使う。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`hprice`のデータを使った上の例を使う。\n",
    "\n",
    "`het_white()`の引数：\n",
    "* 第１引き数：OLSの結果`res_h`の属性`.resid`で誤差項の値\n",
    "* 第２引き数：OLSの結果`res_h`の属性`.model`の属性`exog`を使い定数項を含む説明変数の値"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### 対数変換前"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "het_white(res_h.resid, res_h.model.exog)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "返り値（上から）\n",
    "* `LM statistic`：$LM$統計量\n",
    "* `LM p-value`：$LM$の$p$値\n",
    "* `F statistics`：$F$統計量\n",
    "* `F p-value`：$F$の$p$値\n",
    "\n",
    "一般にWhite検定は$LM$統計量を使ったものを指すが，$F$統計量としても計算できる。\n",
    "\n",
    "5%の有意水準で帰無仮説（$\\text{H}_0$：誤差項は均一分散）を棄却でき，不均一分散の可能性が高い。対処法として変数の変換が挙げられる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### 対数変換後"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "het_white(res_h_log.resid, res_h_log.model.exog)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "5%の有意水準で帰無仮説を棄却できない。即ち，対立仮説（$\\text{H}_A$：帰無仮説は成立しない）を採択し，均一分散の可能性が高い。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "\n",
    "## 残差：図示と線形性"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## 図示"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "仮定４〜６は残差に関するものであり，残差をプロットし不均一分散や非線形性を確認することは回帰分析の重要なステップである。\n",
    "\n",
    "残差を図示する方法として`lmdiag`以外に以下を紹介する。\n",
    "1. `matplotlib`を直接使う\n",
    "2. `seaborn`というパッケージの中にある関数`residplot`を使う\n",
    "\n",
    "上で計算した`res_h`と`res_h_log`を利用し\n",
    "* 横軸：被説明変数の予測値（メソッド`.fittedvalues`）\n",
    "* 縦軸：残差（メソッド`.resid`）\n",
    "\n",
    "となる図を作成する。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### `lmdiag`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "対数変換前"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,7))\n",
    "lmdiag.plot(res_h)\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "対数変換後"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,7))\n",
    "lmdiag.plot(res_h_log)\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "対数変換により残差の変化がより均一的になり，Residuals vs. LeverageのCook's Distanceを見ても外れ値がなくなっている。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### `Matplotlib`の`plot()`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "対数変換前"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "plt.scatter(res_h.fittedvalues, res_h.resid)\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "対数変換後"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "plt.scatter(res_h_log.fittedvalues, res_h_log.resid)\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "対数変換により残差の変化がより均一的になったのが確認できる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### `seaborn`の`residplot()`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`seaborn`は`matplotlib`を利用し様々な図を描ける。`seaborn`については[このサイト](https://seaborn.pydata.org/index.html)を参照。\n",
    "\n",
    "通常`import seaborn as sns`でインポートすることが多いようであるが，ここでは`residplot`のみをインポートしている。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`residplot()`は散布図を作成する関数である。\n",
    "* `x=`：横軸の変数を指定\n",
    "    * 被説明変数を設定することを勧める。\n",
    "* `y=`：縦軸の変数を指定\n",
    "    * 残差\n",
    "* オプション\n",
    "    * `lowerss=True`（デフォルトは`False`）にすると，散布図にベスト・フィットする**曲線**を表示する。\n",
    "\n",
    "---\n",
    "対数変換前"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "residplot(x=res_h.fittedvalues, y=res_h.resid, lowess=True)\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "対数変換後"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "residplot(x=res_h_log.fittedvalues, y=res_h_log.resid, lowess=True)\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## 線形性"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "仮定１で回帰式は線形となっているが，その仮定が正しければ，上の残差の散布図は概ね平行にそして０を中心に上下等間隔に散らばることになる。そのパターンに比較的に近いのは対数変換**後**の図である。不均一分散の場合は，そのパターンが大きく崩れている場合であり，その原因のその１つに回帰式の特定化の間違いがある。例えば，説明変数の２乗が含まれるべきなのに欠落している場合が挙げられる。極端な場合，残差の散布図は$U$字または逆$U$字のような形になりえる。\n",
    "\n",
    "一方，線形性の検定もある。以下ではその１つである RESET (Regression Specification Error Test) 検定を考える。\n",
    "\n",
    "この検定の考え方はそれほど難しくはない。次の回帰式を推定するとしよう。\n",
    "\n",
    "$$y=\\beta_0+\\beta_1x_1+\\beta_2x_2+u$$\n",
    "\n",
    "この線形式が正しければ，$x_1^2$や$x_2^3$等を追加しても統計的に有意ではないはずである。さらに，この線形回帰式の予測値$\\hat{y}$は$x_1$や$x_2$の線形になっているため、$\\hat{y}^2$や$\\hat{y}^3$は$x_1$や$x_2$の非線形となっている。従って，次式を推計し，もし非線形の効果がなければ$\\delta_1$も$\\delta_2$も有意ではないはずである。\n",
    "\n",
    "$$y=\\beta_0+\\beta_1x_1+\\beta_2x_2+\\delta_1\\hat{y}^2+\\delta_2\\hat{y}^3$$\n",
    "\n",
    "この考えに基づいて以下の仮説を検定する。\n",
    "\n",
    "$\\text{H}_0:\\;\\delta_1=\\delta_2=0$（線形回帰式が正しい）\n",
    "\n",
    "$\\text{H}_A$: $\\text{H}_0$は成立しない \n",
    "\n",
    "＜コメント＞\n",
    "* 通常は$y$の3乗まで含めれば十分であろう。\n",
    "* 大標本のもとで$F$検定統計値は$F$分布に従う。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "\n",
    "---\n",
    "`statsmodels`のサブパッケージ`.stats.outliers_influence`のRESET検定用の関数`reset_ramsey`を使う。（`ramsey`はこの検定を考案した学者名前）\n",
    "\n",
    "`reset_ramsey()`の使い方：\n",
    "* 引き数：OLS推定の結果\n",
    "* オプションの`degree`（デフォルトは5）は$y$の何乗までを含めるかを指定する。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "**対数変換前**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "reset_ramsey(res_h,degree=3).summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "返り値\n",
    "* `F`: $F$統計量\n",
    "* `p`: $p$値\n",
    "* `df_denom`: 分母の自由度\n",
    "* `df_num`: 分子の自由度（２乗以上の制約数）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "**対数変換後**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "reset_ramsey(res_h_log,degree=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "5%の有意水準のもとで，`res_h`の帰無仮説を棄却できるが，`res_h_log`では棄却できない。後者の推計式がより適しているようである。"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernel_info": {
   "name": "python3"
  },
  "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.10"
  },
  "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
}