13b_Stochastic_Solow.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 確率的ソロー・モデル"
   ]
  },
  {
   "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 numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.formula.api as smf\n",
    "\n",
    "import japanize_matplotlib\n",
    "import py4macro\n",
    "\n",
    "# numpy v1の表示を使用\n",
    "np.set_printoptions(legacy='1.21')\n",
    "# 警告メッセージを非表示\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## はじめに"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "景気循環を理解する上で全要素生産性の役割が強調される研究があるが，その導入として全要素生産性の特徴について検討する。またソロー・モデルを使い，全要素生産性が確率的に変動する簡単な景気循環モデルのシミュレーションを通してマクロ変数の特徴を考察する。"
   ]
  },
  {
   "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": [
    "マクロ変数の特徴を捉えるにはホワイト・ノイズよりもAR(1)の方がより適している。では，どの様な形でAR(1)をマクロ・モデルに導入すれば良いのだろうか。現在盛んに行われている景気循環研究はDSGEモデル（**D**ynamic **S**tochastic **G**eneral **E**quilibrium Model）と呼ばれるモデルに基づいており、DSGEは名前が示すように次のような特徴がある。\n",
    "* 時間軸に沿ってマクロ変数の変化を考察する動学的な特徴\n",
    "* 確率的な要因の重要な役割\n",
    "* 合理的期待に基づく将来を見据えた消費者と企業の最適行動\n",
    "* 一般均衡のモデル\n",
    "\n",
    "そして，DSGE研究の先駆けとなったモデルが実物的景気循環モデル（**R**eal **B**usiness **C**ycles Model; RBCモデル）と呼ばれる理論である。詳細については後の章に譲るが，基本的なアイデアは簡単で全要素生産性（TFP）がAR(1)に従って確率的に変動し，それが消費者と企業の経済行動に影響することにより景気循環が発生するという理論である。ここではデータを使い全要素生産性の性質を考察する。RBCモデルで需要な役割を果たす労働時間についてもデータを使い特徴を探ることにする。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### データ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "次の生産関数を仮定する。\n",
    "\n",
    "$$\n",
    "Y_t=A_tK_t^aH_t^{1-a}\n",
    "$$\n",
    "\n",
    "* $Y_t$：GDP\n",
    "* $K_t$：資本ストック\n",
    "* $H_t$：総労働時間（労働者数$\\times$労働者平均労働時間）\n",
    "* $A_t$：TFP\n",
    "\n",
    "TFPについて解くと次式を得る。\n",
    "\n",
    "$$\n",
    "A_t=\\dfrac{Y_t}{K_t^aH_t^{1-a}}\n",
    "$$\n",
    "\n",
    "景気循環を考えているため、TFPを計算する際には次の点を考慮することが望ましい。\n",
    "1. 年次データではなく四半期データを使う（年次データでは短期的な変動を捉えることができない）。\n",
    "1. $H_t$は労働者数ではなく，総労働時間とする（短期的な労働供給の変化を捉えるため）。\n",
    "1. $K_t$は資本ストックの稼働率を考慮したデータとする（短期的に実際に生産に使われた資本ストックの変化を捉えるため）。\n",
    "\n",
    "ここでは１と２のみを考慮したデータを使いTFPの変動を考えることにする。これにより資本ストックの稼働率はTFPの一部となってしまうが，RBCモデルの基本的なアイデアは変わらない。\n",
    "\n",
    "データを読み込もう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = py4macro.data('jpn-q')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ここで読み込んだデータ`df`には次の変数が含まれている（詳細は`GDP：水準・トレンド・変動`を参照）。\n",
    "* 期間：1980年Q1〜\n",
    "* `gdp`：国内総生産（支出側）\n",
    "* `capital`：資本ストック\n",
    "* `employed`：就業者\n",
    "* `hours`：労働者1人あたり月平均就業時間（以下では「平均労働時間」と呼ぶ）\n",
    "* `total_hours`：総労働時間（`hours`$\\times$`employed`）\n",
    "\n",
    "まず`gdp`のトレンドからの乖離（％）を`gdp_cycle`として作成しておく。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['gdp_trend_log'] = py4macro.trend( np.log(df['gdp']) )\n",
    "df['gdp_cycle'] = 100 * ( np.log( df['gdp'] ) - df['gdp_trend_log'] )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 労働時間と雇用の特徴"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "平均労働時間`hours`，就業者数`employed`，総労働時間`total_hours`の特徴を考えてみる。まずサイクルを計算するが、次の点に留意し計算方法が異なる。\n",
    "* `hours`は上限があり、長期間持続的に上昇もしくは減少しない。\n",
    "* `employed`は長期間持続的に上昇もしくは減少することは可能な変数。\n",
    "* `total_hours`は、`hours`と`employed`の掛け算となっているため、長期間持続的に上昇もしくは減少することは可能な変数。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "df['hours_trend'] = py4macro.trend(df['hours'])\n",
    "df['hours_cycle'] = 100 * np.log( df['hours']/df['hours_trend'] )\n",
    "\n",
    "df['employed_trend_log'] = py4macro.trend( np.log(df['employed']) )\n",
    "df['employed_cycle'] = 100 * ( np.log( df['employed'] ) - df['employed_trend_log'] )\n",
    "\n",
    "df['total_hours_trend_log'] = py4macro.trend( np.log(df['total_hours']) )\n",
    "df['total_hours_cycle'] = 100 * ( np.log( df['total_hours'] ) - df['total_hours_trend_log'] )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "最初に，平均労働時間を考えよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "df.plot(y=['hours','hours_trend'])\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "平均労働時間は長期的に減少傾向にある。2020年の平均が100に標準化しているので，1980年代には約25％長かったことを示している。次に，就業者数を図示する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "ax_ = df.plot(y='employed')\n",
    "np.exp( df['employed_trend_log'] ).plot(ax=ax_)\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "```{admonition} コードの説明\n",
    "`df['employed_trend_log']`は、トレンドを対数化した値となる。指数化し元の数値に戻すために`np.exp()`を使っている。また、`np.exp(df['employed_trend'])`は`Series`を返すことになり、そのメソッド`plot()`を使って図示している。\n",
    "```\n",
    "\n",
    "一方，雇用者数は増加傾向にある。女性の労働市場参加率の増加や，雇用形態の変化（非正規など）の影響と考えられる。\n",
    "\n",
    "上の２つの変化を反映したのが総労働時間の変化である。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "ax_ = df.plot(y='total_hours')\n",
    "np.exp( df['total_hours_trend_log'] ).plot(ax=ax_)\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "バブル崩壊後は減少傾向にあるが，過去１０年間に持ち直して来ている。これは雇用の拡大の要因が大きい。\n",
    "\n",
    "雇用と労働時間の変動（トレンドからの乖離率）を比べてみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "ax_ = df.plot(y=['hours_cycle','employed_cycle'])\n",
    "ax_.axhline(0, c='red', lw=0.5)\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "次の２つの点を指摘できる。\n",
    "* 労働時間の方が変動が大きい。即ち，総労働時間の調整は就業者数よりも労働時間を通して行われていることが伺える。\n",
    "* 就業者数の変化は、労働時間の変化を後追いする傾向がある。これは就業者の調整費用がより高いためであろう。\n",
    "\n",
    "自己相関係数を計算してみる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "cycle_list = ['hours_cycle','employed_cycle','total_hours_cycle']\n",
    "var_list = ['平均労働時間','就業者数','総労働時間']\n",
    "\n",
    "print('     自己相関係数\\n----------------------')\n",
    "\n",
    "for c, v in zip(cycle_list,var_list):\n",
    "    ac = df[c].autocorr()\n",
    "    print(f'{v:<5}\\t{ac:.3f}')  # 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "```{admonition} コードの説明\n",
    ":class: dropdown\n",
    "\n",
    "1. `:<5`は`f-string`で代入する変数`v`のスペースを`5`に設定し左寄せ`<`にすることを指定している。また`:.3f`は小数点第三位までの表示を指定している。\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "就業者数がより大きな自己相関係数の値となっており，就業者での調整により時間がかかるためである。即ち，今期の就業者数は前期の就業者数に依存するところが大きいという意味である。\n",
    "\n",
    "次に，GDPのトレンドからの乖離との相関係数を計算してみる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "print('GDPサイクルとの相関係数')\n",
    "print('-------------------------')\n",
    "\n",
    "for c, v in zip(cycle_list,var_list):\n",
    "    corr = df[['gdp_cycle',c]].corr().iloc[0,1]\n",
    "    print(f'{v:<9}\\t{corr:.3f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "やはり非常に強い正の相関があり，平均労働時間の影響が強いようである。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[['gdp_cycle', 'total_hours_cycle']].plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 資本ストックの特徴"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['capital_trend'] = py4macro.trend( np.log(df['capital']) )\n",
    "df['capital_cycle'] = 100 * (\n",
    "    np.log( df['capital']) - df['capital_trend'] \n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = np.log(df['capital']).plot()\n",
    "df['capital_trend'].plot(ax=ax)\n",
    "# df['capital_cycle'].plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['capital_cycle'].autocorr()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[['gdp_cycle', 'capital_cycle']].corr().iloc[0,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[['gdp_cycle', 'capital_cycle']].plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### TFPの特徴"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "TFPの水準，トレンド，変動（サイクル）を計算するが，資本の所得分配率を次のように仮定する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "a = 0.36"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# 全要素生産性（対数）の計算\n",
    "df['tfp_log'] = np.log( df['gdp']/( df['capital']**a * df['total_hours']**(1-a) ) )\n",
    "\n",
    "# 全要素生産性のトレンド（対数）の計算\n",
    "df['tfp_log_trend'] = py4macro.trend( df['tfp_log'] )\n",
    "\n",
    "# 全要素生産性のトレンドからの乖離率の計算\n",
    "df['tfp_cycle'] = 100 * ( df['tfp_log'] - df['tfp_log_trend'] )\n",
    "\n",
    "# 全要素生産性とトレンドのプロット\n",
    "ax = df.plot(y=['tfp_log','tfp_log_trend'])\n",
    "ax.set_title('全要素生産性とトレンド（対数）', size=20)\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "全要素生産性の変動をプロットしよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "ax = df.plot(y='tfp_cycle')\n",
    "ax.axhline(0, c='red')\n",
    "ax.set_title('全要素生産性の変動（％）', size=20)\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "この図から全要素生産性の変動は大きいことが分かる。変動の最小値と最大値を確認してみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "df['tfp_cycle'].min(), df['tfp_cycle'].max()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "正の乖離は2％近くあり，負の乖離は4％以上となる。リーマンショック時を除くと正も負も絶対値で概ね2％以内に収まっている。\n",
    "\n",
    "GDPと一緒にプロットしてみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[['gdp_cycle','tfp_cycle']].max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "ax_ = df.plot(y=['gdp_cycle','tfp_cycle'])\n",
    "ax_.set_title('GDPとTFPの変動（％）',size=20)\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "同じ方向に動いているのが確認できる。GDPの変動との相関係数を計算してみると，非常に高いことがわかる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "df[['gdp_cycle','tfp_cycle']].corr().iloc[0,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "以前の計算結果を使うと，この数字はGDPのどの構成要素よりも高いことが確認できる。次に，自己相関係数を計算してみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "df['tfp_cycle'].autocorr()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "`gdp_cycle`程ではないが、強い持続性があることが確認できる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "df['gdp_cycle'].autocorr()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "＜これらの結果の含意＞\n",
    "\n",
    "全要素生産性の変動は大きく持続性も高い。またGDPとの相関性も大きい。これらから全要素生産性が景気循環を引き起こす要因ではないかということが示唆される。全要素生産性は資本や労働時間で説明できないGDPの要素であり，そのような「その他」の部分に景気循環を引き起こすランダムな要素が隠されているかも知れない，という考え方である。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 労働時間とTFPの関係"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "GDPと労働時間には正の相関があり、GDPとTFPにも正の相関がある。\n",
    "では、労働時間とTFPにはどのような関係があるのだろうか。\n",
    "散布図で確認してみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.plot.scatter(x='tfp_cycle', y='total_hours_cycle')\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "正のトレンドがある。相関係数を計算してみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rho_tfp_total_hours = df[['total_hours_cycle', 'tfp_cycle']].corr().iloc[0,1]\n",
    "rho_tfp_total_hours"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "約`0.5`の値となっており、両変数ともGDPと正の相関があるため、驚くべき結果ではないかも知れない。\n",
    "経済学的なメカニズムはどうだろうか？\n",
    "この問は、通常の消費はの労働供給モデルを使うと簡単に理解できる。\n",
    "TFPの上昇は実質賃金の上昇につながるが、次の2つの効果により消費者の労働供給は変化することになる。\n",
    "1. 代替効果：賃金の上昇は余暇の機械費用の上昇となるため、余暇を減らし労働供給を増やす。\n",
    "2. 所得効果：賃金の上昇は所得の増大であり、余暇を増やし労働供給を減らす。\n",
    "\n",
    "労働時間とTFPの正の相関は、代替効果が所得効果を上回っていることを意味する。\n",
    "この点は、TFPを景気循環の要素と考える場合に重要なポイントとなり、[内生的労働供給の導入](sec:14-endogenous_labour)で詳細に議論することになる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### AR(1)としてのTFPと労働時間"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### TFP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "全要素生産性を景気循環の要因として捉えるために，全要素生産性をAR(1)としてモデル化してみようというのがこの節の目的である。そのために`statsmodels`を使い回帰分析をしてみよう。まず１期ずらした変数を作成する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "df['tfp_cycle_lag'] = df['tfp_cycle'].shift()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "```{admonition} コードの説明\n",
    ":class: dropdown\n",
    "\n",
    "`shift()`は`Series`を１期ずらすメソッドである。\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "df[['tfp_cycle','tfp_cycle_lag']].head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "列`tfp_cycle_lag`は列`tfp_cycle`を１期ずらしていることが確認できる。次に回帰分析をおこなう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "res_tfp = smf.ols('tfp_cycle ~ -1 + tfp_cycle_lag', data=df).fit()\n",
    "rho_tfp = res_tfp.params.iloc[0]\n",
    "rho_tfp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "```{admonition} コードの説明\n",
    ":class: dropdown\n",
    "\n",
    "回帰式の中に`-1`とあるが、定数項を省く場合に指定する。\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "自己相関係数と殆ど変わらない値である。２つの変数の散布図を表示してみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "df.plot.scatter(x='tfp_cycle_lag', y='tfp_cycle')\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "正の相関が図でも確認できる。次に$\\epsilon_t$について考える。この変数はホワイト・ノイズであり、TFPショックの「源泉」である。平均0の正規分布としてモデル化するが，上の回帰式の残差の標準偏差をTFPの標準偏差の推定値として考えよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "sigma_tfp = res_tfp.resid.std()\n",
    "sigma_tfp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "```{admonition} コードの説明\n",
    ":class: dropdown\n",
    "\n",
    "* `res_tfp`の属性`.resid`は残差を返している。\n",
    "* `.resid`は`Series`であり，その属性`std()`を使い標本標準偏差を計算している（デフォルトでは`ddof=1`)。\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "初期値は0として`rho_tfp`と`sigma_tfp`を使ってAR(1)の値を計算してプロットしてみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "T = 120        # 標本の大きさ\n",
    "A = 0          # 初期値\n",
    "A_lst = [A]   # 値を格納するリスト\n",
    "\n",
    "for i in range(T):\n",
    "    e = np.random.normal(loc=0, scale=sigma_tfp)\n",
    "    A = rho_tfp * A + e\n",
    "    A_lst.append(A)\n",
    "\n",
    "df_A = pd.Series(A_lst)            # Seriesの作成\n",
    "ax_ = df_A.plot()                   # 図示\n",
    "ax_.axhline(0, c='red')             # 0の平行線\n",
    "\n",
    "ac_tfp = df_A.autocorr()       # 自己相関係数の計算\n",
    "print(f'TFPの自己相関係数：{ac_tfp:.3f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "シミュレーションを行う度に図と自己相関係数は異なることになる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 労働時間"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "労働時間にも同様の分析を行ってみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# `1`期シフトした`total_hours_cycle`を作成\n",
    "df['total_hours_cycle_lag'] = df['total_hours_cycle'].shift()\n",
    "\n",
    "# AR(1)として単回帰分析\n",
    "res_total_hours = smf.ols('total_hours_cycle ~ -1 + total_hours_cycle_lag', data=df).fit()\n",
    "rho_total_hours = res_total_hours.params.iloc[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "スロープ係数の推定値を表示し、散布図をプロットしよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rho_total_hours"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.plot.scatter(x='total_hours_cycle_lag', y='total_hours_cycle')\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "正の相関が確認できる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "source": [
    "#### TFPと労働時間"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "hidden": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "# AR(1)\n",
    "formula1 = 'tfp_cycle ~ -1 + tfp_cycle_lag + total_hours_cycle_lag'\n",
    "formula2 = 'total_hours_cycle ~ -1 + tfp_cycle_lag + total_hours_cycle_lag'\n",
    "res1 = smf.ols(formula1, data=df).fit()\n",
    "res2 = smf.ols(formula2, data=df).fit()\n",
    "print(res1.summary().tables[1])\n",
    "print(res2.summary().tables[1])\n",
    "p1 = res1.params\n",
    "p2 = res2.params\n",
    "np.cov([res1.resid, res2.resid], ddof=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true,
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "# SUR -> the same result\n",
    "from linearmodels.system import SUR\n",
    "df['total_hours_cycle_lag'] = df['total_hours_cycle'].shift()\n",
    "eq1 = 'tfp_cycle ~ tfp_cycle_lag + total_hours_cycle_lag'\n",
    "eq2 = 'total_hours_cycle ~  tfp_cycle_lag + total_hours_cycle_lag'\n",
    "# Create the SUR model\n",
    "model = SUR.from_formula({'eq1': eq1, 'eq2': eq2}, data=df.dropna(subset=['tfp_cycle_lag', 'total_hours_cycle_lag']))\n",
    "results = model.fit()\n",
    "print(results.resids.cov())\n",
    "print(results.summary.tables[1])\n",
    "print(results.summary.tables[2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "#### 相関する2つのAR(1)系列\n",
    "n=16000\n",
    "mean = [0, 0]\n",
    "cov = np.cov([res1.resid, res2.resid], ddof=1)\n",
    "rng = np.random.default_rng()\n",
    "\n",
    "y = 0  # tfp\n",
    "z = 0  # total_hours\n",
    "\n",
    "y_lst = [y]\n",
    "z_lst = [z]\n",
    "uy_lst = [0]\n",
    "uz_lst = [0]\n",
    "\n",
    "for _ in range(n):\n",
    "    u = rng.multivariate_normal(mean, cov)\n",
    "    y = p1['tfp_cycle_lag']*y + p1['total_hours_cycle_lag']*z + u[0]\n",
    "    z = p2['tfp_cycle_lag']*y + p2['total_hours_cycle_lag']*z + u[1]\n",
    "    y_lst.append(y)\n",
    "    z_lst.append(z)\n",
    "    uy_lst.append(u[0])\n",
    "    uz_lst.append(u[1])\n",
    "\n",
    "tmp = pd.DataFrame({'y':y_lst, 'z':z_lst, 'uy':uy_lst, 'uz':uz_lst})\n",
    "\n",
    "tmp['y'].autocorr(), tmp['z'].autocorr(), tmp[['y','z']].corr().iloc[0,1],tmp[['uy','uz']].corr().iloc[0,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 要素分解＝＝＝＝"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma2y = df['gdp_cycle'].var()\n",
    "sigma2a = df['tfp_cycle'].var()\n",
    "sigma2k = df['capital_cycle'].var()\n",
    "sigma2h = df['total_hours_cycle'].var()\n",
    "\n",
    "cov_ak = df[['tfp_cycle', 'capital_cycle']].cov().iloc[0,1]\n",
    "cov_ah = df[['tfp_cycle', 'total_hours_cycle']].cov().iloc[0,1]\n",
    "cov_kh = df[['capital_cycle', 'total_hours_cycle']].cov().iloc[0,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = 0.36\n",
    "\n",
    "# tfp contribution\n",
    "tfp_contribution = ( sigma2a + a*cov_ak + (1-a)*cov_ah ) / sigma2y\n",
    "capital_contribution = ( a**2*sigma2k + a*cov_ak + a*(1-a)*cov_kh ) / sigma2y\n",
    "labour_contribution = ( (1-a)**2*sigma2h + (1-a)*cov_ah + a*(1-a)*cov_kh ) / sigma2y\n",
    "print(tfp_contribution)\n",
    "print(capital_contribution)\n",
    "print(labour_contribution)\n",
    "print(tfp_contribution + capital_contribution + labour_contribution)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = 0.36\n",
    "\n",
    "# tfp contribution\n",
    "tfp_contribution = sigma2a / sigma2y\n",
    "capital_contribution = a**2*sigma2k / sigma2y\n",
    "labour_contribution = (1-a)**2*sigma2h / sigma2y\n",
    "print(tfp_contribution)\n",
    "print(capital_contribution)\n",
    "print(labour_contribution)\n",
    "print(tfp_contribution + capital_contribution + labour_contribution)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# tfp contribution\n",
    "tfp_contribution = sigma2a / sigma2y\n",
    "capital_contribution = sigma2k / sigma2y\n",
    "labour_contribution = sigma2h / sigma2y\n",
    "print(tfp_contribution)\n",
    "print(capital_contribution)\n",
    "print(labour_contribution)\n",
    "print(tfp_contribution + capital_contribution + labour_contribution)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 確率的ソロー・モデル"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### 説明"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "上で日本のデータを使い全要素生産性の`rho`と`sigma`を計算しシミュレーションを行なったが，そのような全要素生産性の変動が景気循環の主な要因と考えるのがRBCモデルと呼ばれるものである。RBCモデルの詳細については後の章で検討するが，ここではその導入として次の点を考察する\n",
    "\n",
    "> AR(1)としてのTFPを組み込んだソロー・モデルを景気循環モデルとして考えシミュレーションをおこない，計算結果をデータと比べてどの程度整合性があるかを検討する。\n",
    "\n",
    "モデルの均衡式として以下を使う。\n",
    "* 生産関数：$Y_t=A_tB_tK_t^aH^{1-a}$\n",
    "    * $A_tB_t$が全要素生産性\n",
    "    * $A_t$は変動を捉える項\n",
    "    \n",
    "        $$\n",
    "        \\begin{align*}\n",
    "        \\log(A_{t})&=\\rho\\log(A_{t-1})+\\varepsilon_t \\\\\n",
    "        \\varepsilon_t&\\sim \\textit{WH}(0,\\sigma^2)\n",
    "        \\end{align*}\n",
    "        $$\n",
    "\n",
    "    * $B_t$はトレンドを捉える項\n",
    "    \n",
    "        $$\n",
    "        B_t=B_0(1+g)^t\n",
    "        $$\n",
    "      \n",
    "    * $H$は一定と仮定する（この仮定によりTFPの効果のみを考察可能となる。）\n",
    "    \n",
    "* 資本の蓄積方程式：$K_{t+1}=sY_t+(1-d)K_t$\n",
    "* 消費：$C_t=(1-s)Y_t$\n",
    "* 投資（貯蓄）：$I_t=sY_t$\n",
    "\n",
    "これらの式から次の2式で$K_t$と$A_t$の値が計算できる。\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "K_{t+1}&=sA_tB_tK_t^aH^{1-a}+(1-d)K_t\\\\\n",
    "\\log(A_{t+1})&=\\rho\\log(A_{t})+\\varepsilon_{t+1}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "変数の初期値を次のように仮定する。\n",
    "* $A_0=1$\n",
    "* $K_0=$ 1980年の資本ストックの平均\n",
    "\n",
    "$H$は次の値として一定とする。\n",
    "* $H=$ 1980年の総労働時間の平均\n",
    "\n",
    "パラメータの値を次を仮定する。\n",
    "* $\\rho$と$\\sigma$は上で計算した`rho_tfp`と`sigma_tfp`の値\n",
    "* $g=$ 1980年から最後の年までの平均四半期成長率\n",
    "* $s=$ 1980年から最後の年までのGDPに占める投資の割合の平均\n",
    "* $d=0.025$：年率約10％を想定\n",
    "* $a=0.36$：資本の所得割合\n",
    "\n",
    "シミュレーションの結果を次の変数を含む`DataFrame`にまとめる。\n",
    "* $Y$，$C$，$I$，$K$，$A$の水準の系列\n",
    "* $Y$，$C$，$I$，$K$，$A$のトレンドの系列\n",
    "* $Y$，$C$，$I$，$K$，$A$の変動（サイクル）の系列"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "### 初期値の計算"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "まずGDP，資本，総労働時間の1980年と2019年の平均を計算する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var_lst = ['gdp', 'capital', 'total_hours']      #1\n",
    "year_lst = [df.index[0].year, df.index[-1].year] #2\n",
    "\n",
    "gdp_dic = {}  #3\n",
    "K_dic = {}    #4\n",
    "H_dic = {}    #5\n",
    "\n",
    "for yr in year_lst:                     #6\n",
    "    \n",
    "    cond = ( df.index.year == yr)        #7\n",
    "    mean = df.loc[cond, var_lst].mean() #8\n",
    "    gdp_dic[yr] = mean.iloc[0]  #9\n",
    "    K_dic[yr] = mean.iloc[1]    #10\n",
    "    H_dic[yr] = mean.iloc[2]    #11"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{admonition} コードの説明\n",
    ":class: dropdown\n",
    "\n",
    "* `#1`：３つの変数のリスト\n",
    "* `#2`：平均を計算する年のリストであり，`df.index[0].year`と`df.index[-1].year`はインデックスから最初と最後の年を抽出する\n",
    "* `#3`：それぞれの年のGDPの平均を格納する辞書\n",
    "* `#4`：それぞれの年の資本ストックの平均を格納する辞書\n",
    "* `#5`：それぞれの年の総労働時間の平均を格納する辞書\n",
    "* `#6`：`year_lst`に対しての`for`ループ\n",
    "* `#7`：`df.index`の年が`yr`と同じ場合は`True`を返し，そうでない場合は`False`を返す`Series`を`cond`に割り当てる。\n",
    "* `#8`：それぞれの年の行，`var_lst`にある列を抽出し平均を計算する。\n",
    "    * `.mean()`は列の平均を計算する。\n",
    "    * 計算した平均は`Series`として返されるので，それを左辺にある変数`mean`に割り当てる。\n",
    "* `#9`：`mean`の0番目の要素を`gdp_dic`のキー`yr`（年）に対応する値をして設定する。\n",
    "* `#10`：`mean`の1番目の要素を`K_dic`のキー`yr`（年）に対応する値をして設定する。\n",
    "* `#11`：`mean`の2番目の要素を`H_dic`のキー`yr`（年）に対応する値をして設定する。\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "このコードの結果を確認してみよう。例えば，"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdp_dic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "from myst_nb import glue\n",
    "var_lst = ['gdp', 'capital', 'total_hours']\n",
    "year_lst = [df.index[0].year, df.index[-1].year]\n",
    "\n",
    "gdp_dic = {}\n",
    "K_dic = {}\n",
    "H_dic = {}\n",
    "\n",
    "for yr in year_lst:\n",
    "    \n",
    "    cond = ( df.index.year == yr)\n",
    "    mean = df.loc[cond, var_lst].mean()\n",
    "    gdp_dic[yr] = mean.iloc[0]\n",
    "    K_dic[yr] = mean.iloc[1]\n",
    "    H_dic[yr] = mean.iloc[2]\n",
    "\n",
    "start_ = list( gdp_dic.keys() )[0]\n",
    "end_ = list( gdp_dic.keys() )[1]\n",
    "glue(\"start_\", start_, display=False)\n",
    "glue(\"end_\", end_, display=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "キー・値の２つのペアがある。キー{glue:}`start_`の値には{glue:}`start_`年のGDPの平均が設定されており，同様に，もう一つのキー{glue:}`end_`の値には{glue:}`end_`年のGDPの平均が設定されている。従って，次のコードで{glue:}`start_`年の値にアクセスできる。\n",
    "```\n",
    "gdp_dic[1980]\n",
    "```\n",
    "他の辞書も同様である。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "````{admonition} なぜforループを使うのか。\n",
    ":class: tip\n",
    "\n",
    "上のコードで`for`ループを使っているが，もし使わなければ次のようなコードになる。\n",
    "```\n",
    "var_lst = ['gdp', 'capital', 'total_hours']\n",
    "mean1980 = df.loc['1980-01-01':'1980-10-01', var_lst].mean()\n",
    "mean2019 = df.loc['2019-01-01':'2019-10-01', var_lst].mean()\n",
    "\n",
    "gdp1980 = mean1980.iloc[0]\n",
    "K1980 = mean1980.iloc[1]\n",
    "H1980 = mean1980.iloc[2]\n",
    "\n",
    "gdp2019 = mean2019.iloc[0]\n",
    "K2019 = mean2019.iloc[1]\n",
    "H2019 = mean2019.iloc[2]\n",
    "```\n",
    "同じようなコードが並んでいる。`py4macro`のデータがアップデートされて`2020`年までのデータが含まれることになったとしよう。`for`ループを使うコードであれば`2019`を`2020`に一ヵ所だけ変更すれば良い（変更する必要もないコードを書くことも可能）。`for`ループを使わないコードで何ヵ所修正する必要があるか数えてみよう。面倒と感じないだろうか。それに修正箇所が増えればbugが紛れ込む可能性が増えていく。\n",
    "````"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "次に全要素生産性の平均四半期成長率を計算するために`df`の列名を確認しよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`gdp`は`0`番目，`capital`は`6`番目，`total_hours`は`11`番目の列にある。この情報をもとに次のように計算しよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var_idx = [0,6,11]         # 1\n",
    "df0 = df.iloc[0,var_idx]   # 2\n",
    "df1 = df.iloc[-1,var_idx]  # 3\n",
    "\n",
    "# 4\n",
    "B0 = df0['gdp'] / ( df0['capital']**a * df0['total_hours']**(1-a) )\n",
    "\n",
    "# 5\n",
    "B1 = df1['gdp'] / ( df1['capital']**a * df1['total_hours']**(1-a) )\n",
    "\n",
    "no_quarters = len(df)      # 6\n",
    "\n",
    "g = ( B1/B0 )**( 1/no_quarters )-1  # 7\n",
    "g"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{admonition} コードの説明\n",
    ":class: dropdown\n",
    "\n",
    "1. 3つの変数の列インデックのリスト\n",
    "2. 最初の行で3つの変数だけを抽出し`df0`に割り当てる。\n",
    "3. 最後の行で3つの変数だけを抽出し`df1`に割り当てる。\n",
    "4. 最初の四半期の全要素生産性の計算\n",
    "5. 最後の四半期の全要素生産性の計算\n",
    "6. 四半期の数\n",
    "7. 平均四半期成長率の計算\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "次にGDPに対する投資の割合の平均を計算する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = ( df['investment'] / df['gdp'] ).mean()\n",
    "s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "約22％が投資に回されている。"
   ]
  },
  {
   "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": [
    "def stochastic_solow(T=160,  # 160=40*4 40年間の四半期の数\n",
    "                     rho=rho_tfp,\n",
    "                     sigma=sigma_tfp,\n",
    "                     g=g,\n",
    "                     s=s,\n",
    "                     d=0.025,\n",
    "                     a=a,\n",
    "                     H=H_dic[1980],  # 1980年平均に固定\n",
    "                     K0=K_dic[1980],\n",
    "                     B0=B0):\n",
    "    \"\"\"引数：\n",
    "            T:     シミュレーションの回数\n",
    "            rho:   AR(1)のパラメータ\n",
    "            sigma: ホワイト・ノイズ分散\n",
    "            g:     TFPの平均四半期成長率\n",
    "            s:     貯蓄率\n",
    "            d:     資本減耗率\n",
    "            a:     生産関数のパラメータ（資本の所得割合）\n",
    "            H:     総労働時間\n",
    "       返り値\n",
    "            次の変数を含むDataFrame\n",
    "            GDP，消費，投資，資本ストック，全要素生産性の水準，トレンド，変動\"\"\"\n",
    "    \n",
    "    # ========== A,B,Kの計算 ==========\n",
    "    # アップデート用変数、最初は初期値を設定\n",
    "    A = 1           #1\n",
    "    B = B0          #2\n",
    "    K = K0          #3\n",
    "    \n",
    "    # 計算結果を格納するリスト\n",
    "    A_lst = [A]     #4\n",
    "    B_lst = [B]     #5\n",
    "    K_lst = [K]     #6\n",
    "\n",
    "    # A,B,Kの時系列の計算\n",
    "    for t in range(1,T):\n",
    "        \n",
    "        K = s * A*B * K**a *H**(1-a) + (1-d)*K  #7\n",
    "        e = np.random.normal(0, scale=sigma)    #8\n",
    "        A = A**rho * np.exp(e)                  #9\n",
    "        B = B0 * (1+g)**t                       #10\n",
    "        \n",
    "        A_lst.append(A)                         #11\n",
    "        B_lst.append(B)\n",
    "        K_lst.append(K)\n",
    "    \n",
    "    # ========== DataFrameの作成 ==========\n",
    "    df_sim = pd.DataFrame({'K':K_lst,           #12\n",
    "                           'A':A_lst,\n",
    "                           'B':B_lst})\n",
    "\n",
    "    # ========== Y,C,Iの計算 ==========          #13\n",
    "    df_sim['Y'] = df_sim['A'] * df_sim['B'] * df_sim['K']**a * H**(1-a)\n",
    "    df_sim['C'] = (1-s) * df_sim['Y']\n",
    "    df_sim['I'] = s * df_sim['Y']\n",
    "    \n",
    "    # ========== トレンドとサイクルの計算 ==========\n",
    "    for v in ['K','A','Y','C','I']:\n",
    "        \n",
    "        # ---------- レンドの計算 ----------      #14\n",
    "        df_sim[f'{v}_log_trend'] = py4macro.trend( np.log(df_sim[v]) )\n",
    "        \n",
    "        # ---------- サイクルの計算 ----------    #15\n",
    "        df_sim[f'{v}の変動'] = 100 * ( np.log( df_sim[v] ) - df_sim[f'{v}_log_trend'] )\n",
    "    \n",
    "    return df_sim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "```{admonition} コードの説明\n",
    ":class: dropdown\n",
    "\n",
    "* `#1`：`for`ループでアップデートされる`A`の値を一次的に格納する。最初は初期値が入っている。\n",
    "* `#2`：`for`ループでアップデートされる`B`の値を一次的に格納する。最初は初期値が入っている。\n",
    "* `#3`：`for`ループでアップデートされる`K`の値を一次的に格納する。最初は初期値が入っている。\n",
    "* `#4`：`A`の値を格納する空のリスト。初期値が入っている。\n",
    "* `#5`：`B`の値を格納する空のリスト。初期値が入っている。\n",
    "* `#6`：`K`の値を格納する空のリスト。初期値が入っている。\n",
    "* `#7`：右辺では(4)~(6)の値（`A`，`B`，`K`）を使い来期の`K`を計算し，その値を左辺の`K`に割り当てる。この瞬間に(6)の`K`がアップデートされる。\n",
    "* `#8`：来期のホワイト・ノイズの値を生成し`e`に割り当てる。\n",
    "* `#9`：右辺では(4)の値（`A`）を使い来期の`A`の値を計算し，その値を左辺の`A`に割り当てる。この瞬間に(4)の`A`がアップデートされる。\n",
    "* `#10`：右辺では来期の`B`の値を計算し，その値を左辺の`B`に割り当てる。この瞬間に(5)の`B`がアップデートされる。\n",
    "* `#11`：リスト(1)~(3)に来期値を追加する。\n",
    "* `#12`：リスト(1)~(3)を使いDataFrameの作成\n",
    "* `#13`：`Y`，`C`，`I`を計算し新たな列としてDataFrameに追加\n",
    "* `#14`：トレンドを計算し新たな列としてDataFrameに追加\n",
    "    * `f-string`を使って新しく作る列名を`K_trend`のように`_trend`を追加する。\n",
    "* `#15`：トレンドからの乖離を計算し新たな列としてDataFrameに追加\n",
    "    * `f-string`を使って新しく作る列名を`K_cycle`のように`_cycle`を追加する。\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### 実行"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "デフォルトの値でシミュレーションを行い、最初の5行を表示する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "df_sim = stochastic_solow()\n",
    "df_sim.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "$Y$の水準とトレンドの図示（対数）。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "np.log( df_sim['Y'] ).plot(title='シミュレーション：Yとトレンド（対数）')\n",
    "df_sim['Y_log_trend'].plot()\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "トレンドの傾きは成長率を表しているが，徐々に減少していることが分かる。最初は成長率に対する資本蓄積の効果が大きいが，定常状態に近づくにつれて資本蓄積の効果が減少しているためである。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "$Y$の変動（サイクル）の図示。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "ax_ = df_sim.plot(y='Yの変動', title='シミュレーション：Yの変動（％）')\n",
    "ax_.axhline(0, c='red')\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "絶対値で概ね2％内に収まった変動となっている。データでは，リーマンショック時はマイナス約4％程乖離したが，それ以外は絶対値で2％以内に収まっており，シミュレーションは概ねデータと同じような特徴を持っていると言える。"
   ]
  },
  {
   "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": [
    "var_lst = ['Y', 'C', 'I', 'K', 'A']\n",
    "\n",
    "print('\\n--- シミュレーション：変動の自己相関係数 ---\\n')\n",
    "\n",
    "for v in var_lst:\n",
    "    \n",
    "    ac = df_sim[f'{v}の変動'].autocorr()\n",
    "    print(f'{v}の変動： {ac:.3f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "まずシミュレーションを行う度にこれらの値は変化することに留意し結果を考えよう。ある程度の自己相関が発生しており，$A$の影響が反映されている。$Y$，$C$，$I$の値は同じなのは，$C$，$I$は$Y$の線形関数であるためであり，これらの３つの変数は$A$の変動と似た動きになっている。また$K$の値が大きいのは、ストック変数であるため変化に時間が掛かることが反映されているためである。\n",
    "\n",
    "シミュレーション結果を評価するためにデータと比べる必要がある。\n",
    "`df`に含まれる日本のデータを使って自己相関係数を計算してみよう。まず消費，投資，資本の変動を計算する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "data_var_lst = ['consumption','investment','capital']\n",
    "\n",
    "for v in data_var_lst:\n",
    "    \n",
    "    df[f'{v}_cycle'] = np.log( df[v]/py4macro.trend(df[v]) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "結果を示すために，`data_var_lst`の最初と最後に`gdp`と`tfp`を追加する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "data_var_lst = ['gdp']+data_var_lst+['tfp']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "print('\\n--- データ：変動の自己相関係数 ---\\n')\n",
    "\n",
    "for v in data_var_lst:\n",
    "    \n",
    "    ac = df[f'{v}_cycle'].autocorr()\n",
    "    print(f'{v:>11}の変動： {ac:>5.3f}')      # 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "```{admonition} コードの説明\n",
    ":class: dropdown\n",
    "\n",
    "1. `f-string`を使って変数`v`と`ac`を代入している。\n",
    "    * `:>11`は`v`を表示する幅を`11`に設定し右寄せ（`>`）にしている。\n",
    "    * `:>5`は`ac`を表示する幅を`5`に設定して右寄せ（`>`）にしている。\n",
    "    * `.3f`は`ac`を小数点第三位までを表示するように指定している。\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "実際の資本ストックは非常に大きな持続性があり、シミュレーションの結果はそれを捉えている。しかし他の変数、特に投資に関しては持続性が足りないようである。\n",
    "\n",
    "しかしシミュレーションの結果はランダム変数の実現値に依存しているので，シミュレーションを行う度に異なる結果が出てくる。また，このシミュレーションでは`T=160`として40年間を考えたが，`T`の値によっても結果は大きく左右される。`T`を大きくするとより安定的な結果になっていくだろう。例えば，`T=10_000`を試してみよう。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true,
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "#### 結果２：GDPとの相関"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "次に、$Y$との相関係数を計算する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "print('\\n--- シミュレーション：GDPとの相関係数 ---\\n')\n",
    "\n",
    "for v in var_lst[1:]:\n",
    "    \n",
    "    cov = df_sim[['Yの変動',f'{v}の変動']].corr().iloc[0,1]\n",
    "    print(f'{v}の変動： {cov:.3f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "* $C$と$I$の値は1.0となっている理由は$Y$と線形になっているためであり，全く同じように変動するためである。\n",
    "* $Y$と$A$は非常に強く相関している。生産関数を見ると，両変数も線形の関係があり，$K$の変化によって変化のズレが生じるメカニズムとなっているためである。\n",
    "\n",
    "シミュレーション結果ををデータと比べてみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "print('\\n--- データ：GDPとの相関係数 ---\\n')\n",
    "\n",
    "for v in data_var_lst[1:]:\n",
    "    cor = df[['gdp_cycle',f'{v}_cycle']].corr().iloc[0,1]\n",
    "    print(f'{v:>11}の変動：{cor:>5.3f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "$C$と$I$に関してはシミュレーションの相関係数は大きすぎる結果となっている。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "#### 総評"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "確率的ソロー・モデルによるシミュレーション結果は，マクロ変数の特徴を捉えている部分もあるが景気循環モデルとして考えるには難しいだろう。消費者の効用最大化問題が欠落しているため，$Y$，$C$，$I$が全く同じように動くのがネックになり，消費や投資の特徴を捉えることができていない。後の章で議論するが，消費者の最適化問題や将来の予測などを導入したRBCモデルは，よりデータに近い結果となる。しかし，いくつか問題が残るのも事実である。\n",
    "* TFPショックとは何なのか？\n",
    "    * 生産関数を見る限りTFPは生産性を表している。技術水準，政府の規制，政治システムや経済制度の変化を捉えると考えられ，例えば，ビジネスに対する規制の強化によって負のショックが発生する。しかし，それだけで実際に観測されるGDPの大きな変動（例えば、リーマン・ショックやコロナ・ショック）を説明するのには疑問である。\n",
    "* ソロー・モデルも含めて価格調整が瞬時に行われるが，次の問題が残る。\n",
    "    * 古典派の二分法が成立すると仮定するが、短期的には無理がある仮定（例えば、賃金）。\n",
    "    * 景気循環に対して金融政策は無力なのか？\n",
    "    \n",
    "このような問題意識を持って，RBCモデルとニューケインジアン・モデルを議論していくことにする。"
   ]
  }
 ],
 "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}