6_Regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(chap:6-regression)=\n",
    "\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 japanize_matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import py4macro\n",
    "import statsmodels.formula.api as smf\n",
    "\n",
    "# numpy v1の表示を使用\n",
    "np.set_printoptions(legacy='1.21')\n",
    "# 警告メッセージを非表示\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "## はじめに"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "この章の目的は２つある。第一に，`Python`を使って回帰分析の方法を紹介する事である。読者には統計学や計量経済学を学び，実際にソフトを使い回帰分析をした経験を持つ人も多く含まれると思う。例えば，無料のRやGretl，そして有料ではあるがStataやEViewsあたりが人気ではないだろうか。`Python`でも回帰分析を簡単におこなうことが可能であり，有名な`statsmodels`と`linearmodels`パッケージを使うと学部の授業で習うことは殆ど全て可能だろう。この章では基本的な回帰分析を考えるので`statsmodels`だけを取り上げ，その基本的なコードを紹介するが，より多くを知りたい場合は[「Pythonで学ぶ入門計量経済学」](https://py4etrics.github.io)を参照して欲しい。\n",
    "\n",
    "この章の第二の目的は，[発展会計](chap:5a-accounting)と[成長会計](chap:5b-accounting)の章で考察した問題を再考することである。その章では，全要素生産性と蓄積生産要素の２つにフォーカスし，それぞれが一人当たりGDPの水準と成長率にどれだけ貢献しているかを考えた。その際，寄与度として分散と共分散を使い計算したが，回帰分析を使うと簡単に計算できる事を示す。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "## `Python`を使った回帰分析 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "[発展会計](chap:5a-accounting)の章で経済間の所得格差の要因を考え，[全要素生産性](sec:5a-tfp)のセクションで全要素生産性と一人当たりGDPの散布図を作成した。そのデータを使い回帰分析のためのコードの書き方を説明する。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true,
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "### データ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "まずデータを作成する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "df2019 = py4macro.data('pwt') \\\n",
    "                 .query('year == 2019') \\\n",
    "                 .reset_index(drop=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "```{admonition} コードの説明\n",
    ":class: dropdown\n",
    "\n",
    "* `py4macro.data()`は`DataFrame`を返すため，そのメソッド`query()`を続けて使っている。また`query()`も`DataFrame`を返すため，そのメソッド`reset_index()`を続けて使っている。\n",
    "* `\\`を使うと改行が可能となる。\n",
    "* `query('year == 2019')`を使って2019年のデータだけを抽出したが，その際に行インデックスが飛び飛びになるので，それを`0`，`1`，`2`....と番号を振り直すために`reset_index()`を使っている。引数`drop=True`がないと，元の飛び飛びの行インデックスが新たな列として残るため，それを残さないようにしている。\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "発展会計の分析で使った同じ変数を作成しよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# 資本の所得シャア\n",
    "a = 1 / 3\n",
    "\n",
    "# 労働者一人当たりGDP\n",
    "df2019['gdp_pc'] = df2019['cgdpo'] / df2019['emp']\n",
    "\n",
    "# 労働者一人当たり資本\n",
    "df2019['k_pc'] = df2019['ck'] / df2019['emp']\n",
    "\n",
    "# 蓄積生産要素\n",
    "df2019['factors'] = df2019['k_pc']**a * ( df2019['avh']*df2019['hc'] )**(1-a)\n",
    "\n",
    "# 全要素生産性\n",
    "df2019['tfp'] = df2019['gdp_pc'] / df2019['factors']\n",
    "\n",
    "# 米国のデータ\n",
    "us2019 = df2019.query('country == \"United States\"')\n",
    "\n",
    "# 労働者一人当たりGDPを標準化（USA=1）\n",
    "df2019['gdp_pc_relative'] = df2019['gdp_pc'] / us2019['gdp_pc'].to_numpy()\n",
    "\n",
    "# 全要素生産性の標準化（USA=1）\n",
    "df2019['tfp_relative'] = df2019['tfp'] / us2019['tfp'].to_numpy()\n",
    "\n",
    "# 蓄積生産要素の標準化（USA=1）\n",
    "df2019['factors_relative'] = df2019['factors'] / us2019['factors'].to_numpy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "全要素生産性と一人当たりGDPの散布図を作成しよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "df2019.plot.scatter(x='tfp_relative', y='gdp_pc_relative')\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 回帰分析"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "上の散布図に直線のトレンドの描く場合\n",
    "\n",
    "$$\n",
    "y_i = a+bA_i\n",
    "$$ (eq:6-regression)\n",
    "\n",
    "となり，`i`は個々の観測値（即ち，経済）を表す。ここでの目的は，散布図のデータを使い定数項`a`とスロープ係数`b`を計算し直線トレンドを描く事である。\n",
    "\n",
    "式[](eq:6-regression)のパラメータ`b`を推計するために`statsmodels`パッケージを使うが，そのサブパッケージ`formula.api`の中に含まれる`ols`を使うことにより、最小二乗法の多くの計算（例えば、検定値など）を自動で行うことが可能である。また推定式を文字列で書くことができるので直感的にコードを書くことが可能となる。より具体的な説明は[このサイト](https://py4etrics.github.io)を参照してほしい。\n",
    "\n",
    "次のステップでコードを書く。\n",
    "1. 文字列で推定式を書く。\n",
    "1. `smf`に含まれる`ols()`関数を使って推定の準備を行う。\n",
    "1. メソッド`.fit()`を使って自動計算を行う。\n",
    "1. 結果の属性やメソッドを使い推定値を表示する。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "**＜ステップ１＞** 推定式の設定\n",
    "\n",
    "回帰式は次のような形で文字列を使って指定する。\n",
    "```\n",
    "'被説明変数 ~ 定数項以外の説明変数'\n",
    "```\n",
    "定数項は自動的に挿入される。また定数項以外の説明変数が複数がる場合は、`+`でつなげるが，今回は単回帰分析となるので説明変数は１つとなる。また回帰式の中で使う変数名は，データが含まれる`DataFrame`の列ラベルを使う。\n",
    "```\n",
    "'gdp_pc_relative ~ tfp_relative'\n",
    "```\n",
    "以下ではこの回帰式を変数`formula`に割り当てているが，変数は分かりやすいものであれば好きなものを使えば良い。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "formula = 'gdp_pc_relative ~ tfp_relative'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "**＜ステップ２＞** 自動計算の準備\n",
    "\n",
    "`smf`の`ols()`関数を使って自動計算の準備として計算の対象となるもの（インスタンスと呼ばれるオブジェクト）を生成する。次のコードでは変数`model`に割り当てる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "model = smf.ols(formula, data=df2019)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`ols()`の第１引数は上で定義した文字列の回帰式であり、第２引数`data`は使用するデータを指定する。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "**＜ステップ３＞** 自動計算\n",
    "\n",
    "`model`のメソッド`.fit()`を使って自動計算し，計算した結果を変数`result`に割り当てる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "result = model.fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "ここでステップ２と３を連続で次のように書いても構わない。\n",
    "```\n",
    "result = ols(formula, data=df).fit()\n",
    "```\n",
    "\n",
    "**＜ステップ４＞** 結果の表示\n",
    "\n",
    "`result`には様々な属性が用意されているが、メソッド`.summary()`を使うと基本的な推定結果を表示することができる。（関数`print()`を使っているが使わなくても同じ情報が表示される。）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "print(result.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "表は3つのセクションから構成されている。\n",
    "* 上段にはOLS推定の基本的な情報が表示されている。\n",
    "    * 左側\n",
    "        * `Dep. Variable`：被説明変数\n",
    "        * `Model`：モデル\n",
    "        * `Method`：手法\n",
    "        * `Data`：日にち\n",
    "        * `Time`：時間\n",
    "        * `No. Observation`：標本の大きさ\n",
    "        * `Df Residuals`：残差の自由度\n",
    "        * `Df Model`：モデルの自由度（定数項以外の説明変数の数）\n",
    "        * `Covariance Type`：共分散のタイプ\n",
    "    * 右側\n",
    "        * `R-squared`：決定係数\n",
    "        * `adj. R-squared`：自由度調整済み決定係数\n",
    "        * `F-statistic`：$F$統計量\n",
    "        * `Prob (F-statistic)`：$F$値\n",
    "        * `Log-Likelihood`：対数尤度\n",
    "        * `AIC`：赤池情報量規準\n",
    "        * `BIC`：ベイズ情報量規準\n",
    "* 中段には主な推定結果が表示される。\n",
    "    * 列ラベル\n",
    "        * `coef`：係数\n",
    "        * `std err`：標準誤差\n",
    "        * `t`：$t$値\n",
    "        * `P>|t|`：$p$値\n",
    "        * `[0.025,0.975]`：信頼区間（5%)\n",
    "    * 行ラベル\n",
    "        * `Intercept`：定数項\n",
    "        * `tfp_relative`：説明変数（選択する変数によって変わる）\n",
    "* 下段には様々な検定などに関する数値が並んでいる。\n",
    "    * 左側\n",
    "        * `Omnibus`：オムニバス検定統計量（帰無仮説：残差は正規分布に従う）\n",
    "        * `Prob(Omnibus)`：オムニバス検定$p$値\n",
    "        * `Skew`：残差の歪度（正規分布であれば`0`）\n",
    "        * `Kurtosis`：残差の尖度（正規分布であれば`3`）\n",
    "    * 右側\n",
    "        * `Durbin-Watson`：ダービン・ワトソン統計量（残差の自己相関の検定）\n",
    "        * `Jarque-Bera (JB)`：ジャーク・ベラ検定統計量（帰無仮説：残差は正規分布に従う）\n",
    "        * `Prob(JB)`：ジャーク・ベラ検定$p$値\n",
    "        * `Cond. No.`：条件指数（Condition Index）の最大値（多重共線性を確認するための指標だが，変数が標準化されて計算していないため変数の値の大きさに依存することになるり，使い難い指標となっている。重回帰分析においての多重共線性の確認については[ここで説明](https://py4etrics.github.io/9_Multiple_Regression.html#sec-9-vif)する手法を使うことを勧める。）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "必要な部分だけを表示する場合は，次のコードを使うと良いだろう。\n",
    "\n",
    "```\n",
    "result.summary().tables[0]\n",
    "result.summary().tables[1]\n",
    "result.summary().tables[2]\n",
    "```\n",
    "\n",
    "係数の推定値に関する中段の表を表示してみる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "print(result.summary().tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "````{note}\n",
    "`statsmodels 13.0`では`.summary()`に引数`slim`（デフォルトは`False`）が追加されており，`.tables`を使わずに次のコードで簡略化された表を表示できる。\n",
    "```\n",
    "result.summary(slim=True)\n",
    "```\n",
    "````"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "結果の解釈を試みてみよう。\n",
    "* 定数項の推定値のp値は小さく、通常の有意水準では「定数項の値はゼロ」の帰無仮説を棄却できる。一方で，一人当たりGDPは正の値を取るが，定数項はマイナスの値となっている。ここでは非線形の関係が示唆されるので，練習問題で再考する。\n",
    "* `tfp_relative`のパラメータの推定値は正の値であり，統計的優位性も高い。推定値を解釈するために，相対的全要素生産性である`tfp_relative`が0.01上昇したとしよう（例えば，0.5から0.51）。相対的一人当たりGDP`gdp_pc_relative`は平均で1.2225x0.01=0.012225増加する事になる。即ち，全要素生産性の１単位の上昇は一人当たりGDPの１単位以上の増加につながっている。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`result`は推定結果に関する情報が詰まったオブジェクトである。何が備わっているか`py4macro`モジュールに含まれる`see()`関数を使って調べてみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "py4macro.see(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "推定結果の表を表示するメソッド`.summary()`が含まれていることが分かる。その他様々なものが含まれているが，詳細は[このサイト](https://www.statsmodels.org/stable/generated/statsmodels.regression.linear_model.RegressionResults.html)（英語）を参考にして欲しい。ここでは代表的なものだけを紹介する。\n",
    "\n",
    "まず係数の推定値は`result`の属性`params`でアクセスできる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "result.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`result.params`のタイプを調べてみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "type(result.params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`Pandas`の`Series`として返されていることが分かる。従って，定数項は`.iloc[]`とインデックス番号を使って`result.params[0]`，スロープ係数は`result.params[1]`で抽出できる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "ahat = result.params.iloc[0]\n",
    "bhat = result.params.iloc[1]\n",
    "ahat, bhat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "もしくは，ラベルを使うこともできる。その場合は，`[]`もしくは`.loc[]`を使うと良い。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "ahat = result.params['Intercept']\n",
    "bhat = result.params['tfp_relative']\n",
    "\n",
    "# これでも同じ結果となる\n",
    "# ahat = result.params.loc['Intercept']\n",
    "# bhat = result.params.loc['tfp_relative']\n",
    "\n",
    "ahat, bhat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "この方法を使い，`tfp_relative`が`x`の場合の`gdp_pc_relative`を計算する関数は作ってみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def calculate_gdp_pc_relative(x):\n",
    "    \n",
    "    gdp = result.params.iloc[0]+result.params.iloc[1]*x\n",
    "    \n",
    "    print(f'相対全要素生産性が{x}の場合の相対的一人当たりGDPは約{gdp:.3f}です。')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`x`が`0.8`の場合を考えよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "calculate_gdp_pc_relative(0.8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "````{note}\n",
    "予測値だけであれば，次のコードでも同じ結果を表示できる。\n",
    "```\n",
    "result.predict({'tfp_relative':0.8})\n",
    "```\n",
    "ここで`.predict()`は予測値を得るメソッド。引数を入れずに使うと，次に説明する属性`.fittedvalues`と同じ値を返す。\n",
    "````"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "次に、標本の散布図に回帰直線を重ねて表示してみる。まず`result`の属性`.fittedvalues`を使い非説明変数の予測値を抽出することができるので、`df2019`に`fitted`のラベルを使って新たな列として追加する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "df2019['OLS fitted'] = result.fittedvalues"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "図を重ねて表示する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "ax_ = df2019.plot.scatter(x='tfp_relative', y='gdp_pc_relative')\n",
    "df2019.sort_values('OLS fitted').plot(x='tfp_relative',\n",
    "                                      y='OLS fitted',\n",
    "                                      color='r', ax=ax_)\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "```{admonition} コードの説明\n",
    ":class: dropdown\n",
    "\n",
    "２行目の`sort_values()`は引数の`'OLS fitted'`を基準に昇順に並べ替える。これを指定しないと直線が綺麗に表示されない。\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "````{note}\n",
    "`gdp_pc_relative`の予測値は，属性`.fittedvalues`もしくはメソッド`.predict()`で取得できるが，次のコードでも同じ結果となる。\n",
    "```\n",
    "df2019['fitted'] = ahat + bhat * df2019.loc[:,'tfp_relative']\n",
    "```\n",
    "````"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "jp-MarkdownHeadingCollapsed": true,
    "tags": [
     "remove-cell"
    ]
   },
   "source": [
    "## 回帰分析による発展会計"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true,
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "(sec:6-regression-da)=\n",
    "\n",
    "## 回帰分析による発展会計"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true,
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "### 説明"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "[発展会計](chap:5a-accounting)の式[](eq:riy-decomposition)では，米国の一人当たりGDPを基準として経済$i$の一人当たりGDPを次式で表した。\n",
    "\n",
    "$$\n",
    "r_i^y =\n",
    "r_i^{\\text{tfp}} +\n",
    "r_i^{\\text{factors}}\n",
    "$$ (eq:6-riy)\n",
    "\n",
    "ここで，\n",
    "* $r_i^y$：一人当たりGDP（対数; 米国＝１）\n",
    "* $r_i^{\\text{tfp}}$：全要素生産性（対数; 米国＝１）\n",
    "* $r_i^{\\text{factors}}$：蓄積生産要素（対数; 米国＝１）\n",
    "\n",
    "また，それぞれの貢献度を次式で数量化した。\n",
    "\n",
    "$$\n",
    "\\text{全要素生産性の寄与度}\n",
    "\\equiv\\beta_{\\text{tfp}}=\n",
    "\\dfrac{\n",
    "    \\text{Var}\\left(r_i^{\\text{tfp}}\\right)\n",
    "    +\\text{Cov}\\left(r_i^{\\text{tfp}},r_i^{\\text{factors}}\\right)\n",
    "    }{\n",
    "    \\text{Var}(r_i^y)\n",
    "    }\n",
    "$$ (eq:6-tfp-contribution)\n",
    "\n",
    "$$\n",
    "\\text{蓄積生産要素の寄与度}\n",
    "\\equiv\n",
    "\\beta_{\\text{factors}}=\n",
    "\\dfrac{\n",
    "    \\text{Var}\\left(r_i^{\\text{factors}}\\right)\n",
    "    +\\text{Cov}\\left(r_i^{\\text{tfp}},r_i^{\\text{factors}}\\right)\n",
    "    }{\n",
    "    \\text{Var}(r_i^y)\n",
    "    }    \n",
    "$$ (eq:6-factors-contribution)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "この章では回帰分析を使って$\\beta_{\\text{tfp}}$と$\\beta_{\\text{factors}}$計算することができることを示す。そのために，まず次の関係が成立することを思い出そう。\n",
    "\n",
    "$$\n",
    "\\text{Cov}(x,x)=\\text{Var}(x)\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\text{Cov}(x,y\\pm z)=\\text{Cov}(x,y)\\pm\\text{Cov}(x,z)\n",
    "$$\n",
    "\n",
    "これらを使うと[](eq:6-tfp-contribution)は次のように書き換えることができる。\n",
    "\n",
    "$$\n",
    "\\beta_{\\text{tfp}}\n",
    "=\\dfrac{\n",
    "    \\text{Var}\\left(r_i^{\\text{tfp}}\\right)\n",
    "    +\\text{Cov}\\left(r_i^{\\text{tfp}},r_i^y-r_i^{\\text{tfp}}\\right)\n",
    "    }{\n",
    "    \\text{Var}(r_i^y)\n",
    "    }\n",
    "=\\dfrac{\n",
    "    \\text{Cov}\\left(r_i^{\\text{tfp}},r_i^{y}\\right)\n",
    "    }{\n",
    "    \\text{Var}(r_i^y)\n",
    "    } \n",
    "$$ (eq:6-tfp-contribution-2)\n",
    "\n",
    "蓄積生産要素の寄与度についても同様に次式で与えられる。\n",
    "\n",
    "$$\n",
    "\\beta_{\\text{factors}}\n",
    "=\\dfrac{\n",
    "    \\text{Cov}\\left(r_i^{\\text{factors}},r_i^{y}\\right)\n",
    "    }{\n",
    "    \\text{Var}(r_i^y)\n",
    "    } \n",
    "$$ (eq:6-factors-contribution-2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "この結果が回帰分析とどのような関係にあるかを確認するために，次の単回帰式を考えよう。\n",
    "\n",
    "$$\n",
    "y_i = a + bx_i+e_i\n",
    "$$\n",
    "\n",
    "最小二乗法を使い推定すると，スロープ係数の推定値は次式で与えられることになる。\n",
    "\n",
    "$$\n",
    "\\hat{b}=\\dfrac{\n",
    "            \\text{Cov}\\left(y,x\\right)\n",
    "            }{\n",
    "            \\text{Var}(x)\n",
    "            } \n",
    "$$ (eq:6-bhat)\n",
    "\n",
    "式[](eq:6-bhat)と２つの式[](eq:6-tfp-contribution-2)と[](eq:6-factors-contribution-2)から，次の回帰式を最小二乗法で推定すると$\\beta_{\\text{tfp}}$と$\\beta_{\\text{factors}}$を計算できることが分かる。\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "&r_i^{\\text{tfp}}=a+\\beta_{\\text{tfp}}r_i^{y}+e_i \\\\\n",
    "&r_i^{\\text{factors}}=c+\\beta_{\\text{factors}}r_i^{y}+u_i\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "以下では`statsmodels`を使って実際に計算してみよう。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 最小二乗法による計算"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`df2019`にある変数を使って，計算に使う変数を作成しよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# 労働者一人当たりGDP（USA=1）の対数化\n",
    "df2019['gdp_pc_relative_log'] = np.log( df2019['gdp_pc_relative'] )\n",
    "\n",
    "# 全要素生産性（USA=1）の対数化\n",
    "df2019['tfp_relative_log'] = np.log( df2019['tfp_relative'] )\n",
    "\n",
    "# 蓄積生産要素（USA=1）の対数化\n",
    "df2019['factors_relative_log'] = np.log( df2019['factors_relative'] )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "次に推定式を定義する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "formula_tfp = 'tfp_relative_log ~ gdp_pc_relative_log'\n",
    "formula_factors = 'factors_relative_log ~ gdp_pc_relative_log'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "自動計算のための準備をし，実際に計算した結果を変数に割り当てる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "result_tfp = smf.ols(formula_tfp, data=df2019).fit()\n",
    "result_factors = smf.ols(formula_factors, data=df2019).fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "パラメータの値を表示してみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "print(f'全要素生産性の寄与度：{result_tfp.params.iloc[1]}\\n'\n",
    "      f'蓄積生産要素の寄与度：{result_factors.params.iloc[1]}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "```{admonition} コードの説明\n",
    ":class: dropdown\n",
    "\n",
    "* `f`は以前説明した`f-string`を使って`{}`に値を代入している。\n",
    "* `()`は複数行にしても構わない。\n",
    "* `''`で囲んだ文字列を複数行にする場合は，それぞれの行で閉じる必要がある。\n",
    "* `\\n`は改行という意味。\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "[全要素生産性と蓄積生産要素の寄与度](sec:5-tfp-factors)で計算した値と同じになることが確認できる。"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}