2018-03-31 Hysteresis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "# 2018-03-31 Hysteresis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting up the Python/Jupyter environment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/javascript": [
       "\n",
       "IPython.OutputArea.prototype._should_scroll = function(lines) {\n",
       "    return false;}"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%javascript\n",
    "\n",
    "IPython.OutputArea.prototype._should_scroll = function(lines) {\n",
    "    return false;}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# keep output cells from shifting to autoscroll: little scrolling\n",
    "# subwindows within the notebook are an annoyance..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/delong/anaconda3/lib/python3.6/site-packages/statsmodels/compat/pandas.py:56: FutureWarning: The pandas.core.datetools module is deprecated and will be removed in a future version. Please use the pandas.tseries module instead.\n",
      "  from pandas.core import datetools\n"
     ]
    }
   ],
   "source": [
    "# set up the environment by reading in every library we might need: \n",
    "# os... graphics... data manipulation... time... math... statistics...\n",
    "\n",
    "import sys\n",
    "import os\n",
    "from urllib.request import urlretrieve\n",
    "\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import Image\n",
    "\n",
    "import pandas as pd\n",
    "from pandas import DataFrame, Series\n",
    "from datetime import datetime\n",
    "\n",
    "import scipy as sp\n",
    "import numpy as np\n",
    "import math\n",
    "import random\n",
    "\n",
    "import seaborn as sns\n",
    "import statsmodels\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "\n",
    "# report library versions..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline \n",
    "\n",
    "# put graphs into the notebook itself..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# graphics setup: seaborn-whitegrid and figure size...\n",
    "\n",
    "plt.style.use('seaborn-whitegrid')\n",
    "\n",
    "figure_size = plt.rcParams[\"figure.figsize\"]\n",
    "figure_size[0] = 12\n",
    "figure_size[1] = 10\n",
    "plt.rcParams[\"figure.figsize\"] = figure_size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# Import series and stuff it in a dictionary\n",
    "\n",
    "# reading in the previously-downloaded long-run real GDP and \n",
    "# GDP per capita file from Johnston and Williamson's \"Measuring \n",
    "# Worth\" website\n",
    "#\n",
    "# constructing a dataframe to hold the data, and then wrapping\n",
    "# that dataframe inside a dict object to append source notes and\n",
    "# links\n",
    "\n",
    "Source_URL = 'http://delong.typepad.com/2018-03-31_q_unemp15-64_realgdp-1.csv' \n",
    "\n",
    "GDPUN1564_df = pd.read_csv(\n",
    "    Source_URL, \n",
    "    parse_dates = True,\n",
    "    index_col = 0)\n",
    "\n",
    "GDPUN1564_dict = {}\n",
    "GDPUN1564_dict[\"Source\"] = \"FRED: http://fred.stlouisfed.org\"\n",
    "GDPUN1564_dict[\"SourceURL\"] = Source_URL\n",
    "\n",
    "GDPUN1564_dict[\"SourceNotes\"] = \"Quarterly: Real GDP (2009) and 15-64 Unemployment Rate\"\n",
    "\n",
    "GDPUN1564_dict[\"SourceDescription\"] = \"https://fred.stlouisfed.org/graph/?id=LRUN64TTUSQ156S, https://fred.stlouisfed.org/graph/?id=GDPC1,\"\n",
    "\n",
    "\n",
    "GDPUN1564_dict[\"df\"] = GDPUN1564_df\n",
    "GDPUN1564_df.head()\n",
    "\n",
    "# create an observation number variable and a\n",
    "# log real GDP variable in the GDPUN1564_df\n",
    "\n",
    "GDPUN1564_df.insert(0, 'Observation', range(0, len(GDPUN1564_df)))\n",
    "GDPUN1564_df[\"LN_RGDP\"] = np.log(GDPUN1564_df.RGDP)\n",
    "\n",
    "GDPUN1564_df[\"UN_lag\"] = GDPUN1564_df.UNEMP1564.shift(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:              UNEMP1564   R-squared:                       0.564\n",
      "Model:                            OLS   Adj. R-squared:                  0.562\n",
      "Method:                 Least Squares   F-statistic:                     241.0\n",
      "Date:                Fri, 06 Apr 2018   Prob (F-statistic):           2.12e-35\n",
      "Time:                        15:13:55   Log-Likelihood:                -273.57\n",
      "No. Observations:                 188   AIC:                             551.1\n",
      "Df Residuals:                     186   BIC:                             557.6\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          1.5478      0.323      4.796      0.000       0.911       2.185\n",
      "UN_lag         0.7572      0.049     15.523      0.000       0.661       0.853\n",
      "==============================================================================\n",
      "Omnibus:                       60.884   Durbin-Watson:                   0.157\n",
      "Prob(Omnibus):                  0.000   Jarque-Bera (JB):              119.856\n",
      "Skew:                           1.577   Prob(JB):                     9.41e-27\n",
      "Kurtosis:                       5.315   Cond. No.                         28.7\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "\n",
    "Unemployment = GDPUN1564_df.UNEMP1564[4:]\n",
    "X = GDPUN1564_df.UN_lag[4:]\n",
    "Unemployment_Lagged_1_Year = sm.add_constant(X)\n",
    "results = smf.OLS(Unemployment, Unemployment_Lagged_1_Year).fit()\n",
    "print(results.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:              UNEMP1564   R-squared:                       0.552\n",
      "Model:                            OLS   Adj. R-squared:                  0.549\n",
      "Method:                 Least Squares   F-statistic:                     177.3\n",
      "Date:                Fri, 06 Apr 2018   Prob (F-statistic):           7.08e-27\n",
      "Time:                        15:14:52   Log-Likelihood:                -199.57\n",
      "No. Observations:                 146   AIC:                             403.1\n",
      "Df Residuals:                     144   BIC:                             409.1\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          1.5707      0.361      4.352      0.000       0.857       2.284\n",
      "UN_lag         0.7476      0.056     13.317      0.000       0.637       0.859\n",
      "==============================================================================\n",
      "Omnibus:                       43.139   Durbin-Watson:                   0.186\n",
      "Prob(Omnibus):                  0.000   Jarque-Bera (JB):               73.594\n",
      "Skew:                           1.445   Prob(JB):                     1.05e-16\n",
      "Kurtosis:                       4.934   Cond. No.                         30.0\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "\n",
    "Unemployment = GDPUN1564_df.UNEMP1564[4:150]\n",
    "X = GDPUN1564_df.UN_lag[4:150]\n",
    "Unemployment_Lagged_1_Year = sm.add_constant(X)\n",
    "results = smf.OLS(Unemployment, Unemployment_Lagged_1_Year).fit()\n",
    "print(results.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:              UNEMP1564   R-squared:                       0.535\n",
      "Model:                            OLS   Adj. R-squared:                  0.523\n",
      "Method:                 Least Squares   F-statistic:                     44.86\n",
      "Date:                Fri, 06 Apr 2018   Prob (F-statistic):           5.58e-08\n",
      "Time:                        15:14:46   Log-Likelihood:                -68.941\n",
      "No. Observations:                  41   AIC:                             141.9\n",
      "Df Residuals:                      39   BIC:                             145.3\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          1.8016      0.806      2.235      0.031       0.171       3.432\n",
      "UN_lag         0.7413      0.111      6.698      0.000       0.517       0.965\n",
      "==============================================================================\n",
      "Omnibus:                       15.916   Durbin-Watson:                   0.097\n",
      "Prob(Omnibus):                  0.000   Jarque-Bera (JB):               17.727\n",
      "Skew:                           1.478   Prob(JB):                     0.000141\n",
      "Kurtosis:                       4.278   Cond. No.                         28.7\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "\n",
    "Unemployment = GDPUN1564_df.UNEMP1564[151:]\n",
    "X = GDPUN1564_df.UN_lag[151:]\n",
    "Unemployment_Lagged_1_Year = sm.add_constant(X)\n",
    "results = smf.OLS(Unemployment, Unemployment_Lagged_1_Year).fit()\n",
    "print(results.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>        <td>UNEMP1564</td>    <th>  R-squared:         </th> <td>   0.535</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.523</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   44.86</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Fri, 06 Apr 2018</td> <th>  Prob (F-statistic):</th> <td>5.58e-08</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>13:48:29</td>     <th>  Log-Likelihood:    </th> <td> -68.941</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>    41</td>      <th>  AIC:               </th> <td>   141.9</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>    39</td>      <th>  BIC:               </th> <td>   145.3</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     1</td>      <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "     <td></td>       <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th>  <td>    1.8016</td> <td>    0.806</td> <td>    2.235</td> <td> 0.031</td> <td>    0.171</td> <td>    3.432</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>UN_lag</th> <td>    0.7413</td> <td>    0.111</td> <td>    6.698</td> <td> 0.000</td> <td>    0.517</td> <td>    0.965</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>15.916</td> <th>  Durbin-Watson:     </th> <td>   0.097</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.000</td> <th>  Jarque-Bera (JB):  </th> <td>  17.727</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td> 1.478</td> <th>  Prob(JB):          </th> <td>0.000141</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 4.278</td> <th>  Cond. No.          </th> <td>    28.7</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:              UNEMP1564   R-squared:                       0.535\n",
       "Model:                            OLS   Adj. R-squared:                  0.523\n",
       "Method:                 Least Squares   F-statistic:                     44.86\n",
       "Date:                Fri, 06 Apr 2018   Prob (F-statistic):           5.58e-08\n",
       "Time:                        13:48:29   Log-Likelihood:                -68.941\n",
       "No. Observations:                  41   AIC:                             141.9\n",
       "Df Residuals:                      39   BIC:                             145.3\n",
       "Df Model:                           1                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const          1.8016      0.806      2.235      0.031       0.171       3.432\n",
       "UN_lag         0.7413      0.111      6.698      0.000       0.517       0.965\n",
       "==============================================================================\n",
       "Omnibus:                       15.916   Durbin-Watson:                   0.097\n",
       "Prob(Omnibus):                  0.000   Jarque-Bera (JB):               17.727\n",
       "Skew:                           1.478   Prob(JB):                     0.000141\n",
       "Kurtosis:                       4.278   Cond. No.                         28.7\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "\n",
    "Unemp = GDPUN1564_df.UNEMP1564[151:]\n",
    "Un_lag = GDPUN1564_df.UN_lag[151:]\n",
    "Un_lag = sm.add_constant(Un_lag)\n",
    "results = smf.OLS(Unemp, Un_lag).fit()\n",
    "results.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# convergence to the balanced growth path\n",
    "#\n",
    "# we need to alter our dataframe in order to add a BGP line\n",
    "#\n",
    "# # we are going to want to see what happens for lots of\n",
    "# different model parameter values and initial conditions,\n",
    "# so stuff our small simulation program inside a function, so \n",
    "# we can then invoke it with a single line...\n",
    "\n",
    "def sgm_bgp_100yr_run(L0, E0, K0, n=0.01, \n",
    "    g=0.02, s=0.15, alpha=0.5, delta=0.03, \n",
    "    T = 100):\n",
    "\n",
    "    sg_df = pd.DataFrame(index=range(T),columns=['Labor', \n",
    "        'Efficiency',\n",
    "        'Capital',\n",
    "        'Output',\n",
    "        'Output_per_Worker',\n",
    "        'Capital_Output_Ratio',\n",
    "        'BGP_Output',\n",
    "        'BGP_Output_per_Worker',\n",
    "        'BGP_Capital_Output_Ratio',\n",
    "        'BGP_Capital'],\n",
    "        dtype='float')\n",
    "\n",
    "    sg_df.Labor[0] = L0\n",
    "    sg_df.Efficiency[0] = E0\n",
    "    sg_df.Capital[0] = K0\n",
    "    sg_df.Output[0] = (sg_df.Capital[0]**alpha * (sg_df.Labor[0] * \n",
    "        sg_df.Efficiency[0])**(1-alpha))\n",
    "    sg_df.Output_per_Worker[0] = sg_df.Output[0]/sg_df.Labor[0]\n",
    "    sg_df.Capital_Output_Ratio[0] = sg_df.Capital[0]/sg_df.Output[0]\n",
    "    \n",
    "    sg_df.BGP_Capital_Output_Ratio[0] = (s / (n + g + delta))\n",
    "    sg_df.BGP_Output_per_Worker[0] = sg_df.Efficiency[0] * (\n",
    "        sg_df.BGP_Capital_Output_Ratio[0]*(alpha/(1 - alpha)))\n",
    "    sg_df.BGP_Output[0] = sg_df.BGP_Output_per_Worker[0] * sg_df.Labor[0]\n",
    "    sg_df.BGP_Capital[0] = sg_df.Labor[0] * sg_df.Efficiency[0] * (\n",
    "        sg_df.BGP_Capital_Output_Ratio[0]*(1/(1 - alpha)))\n",
    "\n",
    "    for i in range(T):\n",
    "        sg_df.Labor[i+1] = sg_df.Labor[i] + sg_df.Labor[i] * n\n",
    "        sg_df.Efficiency[i+1] = sg_df.Efficiency[i] + sg_df.Efficiency[i] * g\n",
    "        sg_df.Capital[i+1] = sg_df.Capital[i] - sg_df.Capital[i] * delta + (\n",
    "            sg_df.Output[i] * s)\n",
    "        sg_df.Output[i+1] = (sg_df.Capital[i+1]**alpha * \n",
    "            (sg_df.Labor[i+1] * sg_df.Efficiency[i+1])**(1-alpha))\n",
    "        sg_df.Output_per_Worker[i+1] = sg_df.Output[i+1]/sg_df.Labor[i+1]\n",
    "        sg_df.Capital_Output_Ratio[i+1] = (sg_df.Capital[i+1]/\n",
    "            sg_df.Output[i+1])\n",
    "        sg_df.BGP_Capital_Output_Ratio[i+1] = (s / (n + g + delta))\n",
    "        sg_df.BGP_Output_per_Worker[i+1] = sg_df.Efficiency[i+1] * (\n",
    "            sg_df.BGP_Capital_Output_Ratio[i+1]**(alpha/(1 - alpha)))\n",
    "        sg_df.BGP_Output[i+1] = (sg_df.BGP_Output_per_Worker[i+1] * \n",
    "            sg_df.Labor[i+1])\n",
    "        sg_df.BGP_Capital[i+1] = (s / (n + g + delta))**(1/(1-alpha)) * (\n",
    "            sg_df.Efficiency[i+1] * sg_df.Labor[i+1])\n",
    "        \n",
    "    fig = plt.figure(figsize=(12, 12))\n",
    "\n",
    "    ax1 = plt.subplot(3,2,1)\n",
    "    sg_df.Labor.plot(ax = ax1, title = \"Labor Force\")\n",
    "    plt.ylabel(\"Parameters\")\n",
    "    plt.ylim(0, )\n",
    "\n",
    "    ax2 = plt.subplot(3,2,2)\n",
    "    sg_df.Efficiency.plot(ax = ax2, title = \"Efficiency of Labor\")\n",
    "    plt.ylim(0, )\n",
    "    \n",
    "    ax3 = plt.subplot(3,2,3)\n",
    "    sg_df.BGP_Capital.plot(ax = ax3, title = \"BGP Capital Stock\")\n",
    "    sg_df.Capital.plot(ax = ax3, title = \"Capital Stock\")\n",
    "    plt.ylabel(\"Values\")\n",
    "    plt.ylim(0, )\n",
    "\n",
    "    ax4 = plt.subplot(3,2,4)\n",
    "    sg_df.BGP_Output.plot(ax = ax4, title = \"BGP Output\")\n",
    "    sg_df.Output.plot(ax = ax4, title = \"Output\")\n",
    "    plt.ylim(0, )\n",
    "\n",
    "    ax5 = plt.subplot(3,2,5)\n",
    "    sg_df.BGP_Output_per_Worker.plot(ax = ax5, \n",
    "        title = \"BGP Output per Worker\")\n",
    "    sg_df.Output_per_Worker.plot(ax = ax5, title = \"Output per Worker\")\n",
    "    plt.xlabel(\"Years\")\n",
    "    plt.ylabel(\"Ratios\")\n",
    "    plt.ylim(0, )\n",
    "\n",
    "    ax6 = plt.subplot(3,2,6)\n",
    "    sg_df.BGP_Capital_Output_Ratio.plot(ax = ax6, \n",
    "        title = \"BGP Capital-Output Ratio\")\n",
    "    sg_df.Capital_Output_Ratio.plot(ax = ax6, \n",
    "        title = \"Capital-Output Ratio\")\n",
    "    plt.xlabel(\"Years\")\n",
    "    plt.ylim(0, )\n",
    "\n",
    "    plt.suptitle('Solow Growth Model: Simulation Run', size = 20)\n",
    "\n",
    "    plt.show()\n",
    "    \n",
    "    print(n, \"is the labor force growth rate\")\n",
    "    print(g, \"is the efficiency of labor growth rate\")\n",
    "    print(delta, \"is the depreciation rate\")\n",
    "    print(s, \"is the savings rate\")\n",
    "    print(alpha, \"is the decreasing-returns-to-scale parameter\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "sgm_bgp_100yr_run(1000, 1, 100, n=0.005, g=0.01, \n",
    "    s=0.225, alpha=0.5, delta=0.03)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# convergence to the balanced growth path—log graphs\n",
    "#\n",
    "# we need to alter our dataframe in order to add a BGP line\n",
    "#\n",
    "# # we are going to want to see what happens for lots of\n",
    "# different model parameter values and initial conditions,\n",
    "# so stuff our small simulation program inside a function, so \n",
    "# we can then invoke it with a single line...\n",
    "\n",
    "def log_sgm_bgp_100yr_run(L0, E0, K0, n=0.01, g=0.02, \n",
    "    s=0.15, alpha=0.5, delta=0.03, T=100):\n",
    "\n",
    "    sg_df = pd.DataFrame(index=range(T),columns=['Labor', \n",
    "        'Efficiency',\n",
    "        'Capital',\n",
    "        'Output',\n",
    "        'Output_per_Worker',\n",
    "        'Capital_Output_Ratio',\n",
    "        'BGP_Output',\n",
    "        'BGP_Output_per_Worker',\n",
    "        'BGP_Capital_Output_Ratio',\n",
    "        'BGP_Capital'],\n",
    "        dtype='float')\n",
    "\n",
    "    sg_df.Labor[0] = L0\n",
    "    sg_df.Efficiency[0] = E0\n",
    "    sg_df.Capital[0] = K0\n",
    "    sg_df.Output[0] = (sg_df.Capital[0]**alpha * (sg_df.Labor[0] * sg_df.Efficiency[0])**(1-alpha))\n",
    "    sg_df.Output_per_Worker[0] = sg_df.Output[0]/sg_df.Labor[0]\n",
    "    sg_df.Capital_Output_Ratio[0] = sg_df.Capital[0]/sg_df.Output[0]\n",
    "    sg_df.BGP_Capital_Output_Ratio[0] = (s / (n + g + delta))\n",
    "    sg_df.BGP_Output_per_Worker[0] = sg_df.Efficiency[0] * sg_df.BGP_Capital_Output_Ratio[0]*(alpha/(1 - alpha))\n",
    "    sg_df.BGP_Output[0] = sg_df.BGP_Output_per_Worker[0] * sg_df.Labor[0]\n",
    "    sg_df.BGP_Capital[0] = (s / (n + g + delta))**(1/(1-alpha)) * sg_df.Efficiency[0] * sg_df.Labor[0]\n",
    "\n",
    "    for i in range(T):\n",
    "        sg_df.Labor[i+1] = sg_df.Labor[i] + sg_df.Labor[i] * n\n",
    "        sg_df.Efficiency[i+1] = sg_df.Efficiency[i] + sg_df.Efficiency[i] * g\n",
    "        sg_df.Capital[i+1] = sg_df.Capital[i] - sg_df.Capital[i] * delta + sg_df.Output[i] * s \n",
    "        sg_df.Output[i+1] = (sg_df.Capital[i+1]**alpha * (sg_df.Labor[i+1] * sg_df.Efficiency[i+1])**(1-alpha))\n",
    "        sg_df.Output_per_Worker[i+1] = sg_df.Output[i+1]/sg_df.Labor[i+1]\n",
    "        sg_df.Capital_Output_Ratio[i+1] = sg_df.Capital[i+1]/sg_df.Output[i+1]\n",
    "        sg_df.BGP_Capital_Output_Ratio[i+1] = (s / (n + g + delta))\n",
    "        sg_df.BGP_Output_per_Worker[i+1] = sg_df.Efficiency[i+1] * sg_df.BGP_Capital_Output_Ratio[i+1]**(alpha/(1 - alpha))\n",
    "        sg_df.BGP_Output[i+1] = sg_df.BGP_Output_per_Worker[i+1] * sg_df.Labor[i+1]\n",
    "        sg_df.BGP_Capital[i+1] = (s / (n + g + delta))**(1/(1-alpha)) * sg_df.Efficiency[i+1] * sg_df.Labor[i+1]\n",
    "\n",
    "        \n",
    "    fig = plt.figure(figsize=(12, 12))\n",
    "\n",
    "    ax1 = plt.subplot(3,2,1)\n",
    "    np.log(sg_df.Labor).plot(ax = ax1, title = \"Labor Force\")\n",
    "    plt.ylabel(\"Log Values\")\n",
    "    plt.ylim(0, )\n",
    "\n",
    "    ax2 = plt.subplot(3,2,2)\n",
    "    np.log(sg_df.Efficiency).plot(ax = ax2, title = \"Efficiency of Labor\")\n",
    "    plt.ylim(0, )\n",
    "    \n",
    "    ax3 = plt.subplot(3,2,3)\n",
    "    np.log(sg_df.BGP_Capital).plot(ax = ax3, title = \"BGP Capital Stock\")\n",
    "    np.log(sg_df.Capital).plot(ax = ax3, title = \"Capital Stock\")\n",
    "    plt.ylabel(\"Log Values\")\n",
    "    plt.ylim(0, )\n",
    "\n",
    "    ax4 = plt.subplot(3,2,4)\n",
    "    np.log(sg_df.BGP_Output).plot(ax = ax4, title = \"BGP Output\")\n",
    "    np.log(sg_df.Output).plot(ax = ax4, title = \"Output\")\n",
    "    plt.ylim(0, )\n",
    "\n",
    "    ax5 = plt.subplot(3,2,5)\n",
    "    np.log(sg_df.BGP_Output_per_Worker).plot(ax = ax5, title = \"BGP Output per Worker\")\n",
    "    np.log(sg_df.Output_per_Worker).plot(ax = ax5, title = \"Output per Worker\")\n",
    "    plt.xlabel(\"Years\")\n",
    "    plt.ylabel(\"Log Ratios\")\n",
    "    plt.ylim(0, )\n",
    "\n",
    "    ax6 = plt.subplot(3,2,6)\n",
    "    np.log(sg_df.BGP_Capital_Output_Ratio).plot(ax = ax6, title = \"BGP Capital-Output Ratio\")\n",
    "    np.log(sg_df.Capital_Output_Ratio).plot(ax = ax6, title = \"Capital-Output Ratio\")\n",
    "    plt.xlabel(\"Years\")\n",
    "    plt.ylim(0, )\n",
    "\n",
    "    plt.suptitle('Solow Growth Model: Simulation Run: Log Scale', size = 20)\n",
    "\n",
    "    plt.show()\n",
    "    \n",
    "    print(n, \"is the labor force growth rate\")\n",
    "    print(g, \"is the efficiency of labor growth rate\")\n",
    "    print(delta, \"is the depreciation rate\")\n",
    "    print(s, \"is the savings rate\")\n",
    "    print(alpha, \"is the decreasing-returns-to-scale parameter\")\n",
    "    \n",
    "log_sgm_bgp_100yr_run(1000, 1, 100, n=0.005, g=0.01, s=0.24, \n",
    "    alpha=0.5, delta=0.03)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# suppose we started the economy on some balanced growth path, say\n",
    "# for s = 0.20. And then s jumped to 0.25. What would happen?\n",
    "#\n",
    "# n=0.01, g=0.01, delta=0.03, s=0.20, alpha=0.5...\n",
    "# SS K/Y = 4...\n",
    "# Y/L = 4 x E\n",
    "# K/L = 16 x E\n",
    "#\n",
    "# start the economy on its balanced growth path..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "log_sgm_bgp_100yr_run(1000, 1, 16000, n=0.01, g=0.01, s=0.20, \n",
    "    alpha=0.5, delta=0.03)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# yup...\n",
    "# we can look at it in levels too:\n",
    "\n",
    "sgm_bgp_100yr_run(1000, 1, 16000, n=0.01, g=0.01, \n",
    "    s=0.20, alpha=0.5, delta=0.03)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# Now, from the s = 0.20 BGP, what happens if we suddenly jump s to 0.25\n",
    "# and keep it there?\n",
    "#\n",
    "# This happens:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# in levels:\n",
    "\n",
    "sgm_bgp_100yr_run(1000, 1, 16000, n=0.01, g=0.01, \n",
    "    s=0.25, alpha=0.5, delta=0.03)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "log_sgm_bgp_100yr_run(1000, 1, 16000, n=0.01, g=0.01, s=0.25, \n",
    "    alpha=0.5, delta=0.03)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# in the U.S. today:\n",
    "\n",
    "n = 0.01\n",
    "g = 0.015\n",
    "delta = 0.03\n",
    "\n",
    "sinitial = 0.22\n",
    "alpha = 0.333\n",
    " \n",
    "KoverYstarinitial = sinitial/(n + g + delta)\n",
    "\n",
    "E = 65068\n",
    "\n",
    "YoverLstarinitial = KoverYstarinitial**(alpha/(1-alpha)) * E\n",
    "\n",
    "print(KoverYstarinitial, \"= KoverYstarinitial\")\n",
    "print(YoverLstarinitial, \"= YoverLstarinitial\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# a \"what if\"—if the tax \"reform\" were to boost the savings rate by\n",
    "# 1.4% points...\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "n = 0.01\n",
    "g = 0.015\n",
    "delta = 0.03\n",
    "\n",
    "sfinal = 0.234\n",
    "alpha = 0.333\n",
    " \n",
    "KoverYstarfinal = sfinal/(n + g + delta)\n",
    "\n",
    "E = 65068\n",
    "\n",
    "YoverLstarfinal = KoverYstarfinal**(alpha/(1-alpha)) * E\n",
    "long_run_growth_effect = np.log(YoverLstarfinal/YoverLstarinitial)\n",
    "\n",
    "print(KoverYstarfinal, \"= KoverYstarfinal\")\n",
    "print(YoverLstarfinal, \"= YoverLstarfinal\")\n",
    "print(np.log(YoverLstarfinal/YoverLstarinitial), \"= long-run growth effect\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# speed of convergence\n",
    "\n",
    "speed_of_convergence = ((1 - alpha)*\n",
    "    (n+g+delta))\n",
    "\n",
    "print(speed_of_convergence, \n",
    "    \"= the speed of convergence\")\n",
    "\n",
    "initial_year_growth_boost = (long_run_growth_effect \n",
    "    * speed_of_convergence)\n",
    "\n",
    "print(initial_year_growth_boost, \n",
    "    \"= initial year growth boost\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "**Quoting the four Stanford economists (plus five others):**\n",
    "\n",
    "<img style=\"display:block; margin-left:auto; margin-right:auto;\" src=\"http://delong.typepad.com/.a/6a00e551f08003883401bb09ef5382970d-pi\" alt=\"A conventional approach\" title=\"a_conventional_approach.png\" border=\"0\" width=\"599\" height=\"197\" />\n",
    "\n",
    "<https://www.wsj.com/articles/how-tax-reform-will-lift-the-economy-1511729894?mg=prod/accounts-wsj>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# what is (1 - alpha)(n + g + delta) here?\n",
    "#\n",
    "# (1 - alpha) = 2/3\n",
    "# (n + g + delta) = .045\n",
    "\n",
    "convergence_speed = 2/3 * .045\n",
    "\n",
    "print(convergence_speed, \"= convergence speed\")\n",
    "\n",
    "for i in range(11):\n",
    "    print(.03 - .03 * np.exp(-convergence_speed * i), \"= growth over\", i, \"years\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "##### what if we raise alpha to the \"DeLong Summers\" value?\n",
    "#\n",
    "# note: Brad DeLong and Larry Summers do **not** believe\n",
    "# that the tax \"reform\" will raise the savings rate by \n",
    "# 1.34%: this is a \"what if\" exercise that, as far as I\n",
    "# know, is not the analytical position of anybody:\n",
    "\n",
    "n = 0.01\n",
    "g = 0.015\n",
    "delta = 0.03\n",
    "\n",
    "sinitial = 0.22\n",
    "sfinal = 0.234\n",
    "alpha = 0.55\n",
    "\n",
    "KoverYstarinitial = sinitial/(n + g + delta)\n",
    "KoverYstarfinal = sfinal/(n + g + delta)\n",
    "\n",
    "E = 22000\n",
    "\n",
    "YoverLstarfinal = KoverYstarfinal**(alpha/(1-alpha)) * E\n",
    "YoverLstarinitial = KoverYstarinitial**(alpha/(1-alpha)) * E\n",
    "\n",
    "print(KoverYstarfinal, \"= KoverYstarfinal\")\n",
    "print(YoverLstarfinal, \"= YoverLstarfinal\")\n",
    "print(np.log(YoverLstarfinal/YoverLstarinitial), \n",
    "    \"= long-run growth effect\")\n",
    "print((1 - alpha)*(n+g+delta), \"= speed of convergence\")\n",
    "print((1 - alpha)*(n+g+delta)*\n",
    "    np.log(YoverLstarfinal/YoverLstarinitial),\n",
    "    \"= first year growth effect\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# raising the alpha parameter even further...\n",
    "\n",
    "n = 0.01\n",
    "g = 0.015\n",
    "delta = 0.03\n",
    "\n",
    "sinitial = 0.22\n",
    "sfinal = 0.234\n",
    "alpha = 0.667\n",
    "\n",
    "KoverYstarinitial = sinitial/(n + g + delta)\n",
    "KoverYstarfinal = sfinal/(n + g + delta)\n",
    "\n",
    "E = 7151\n",
    "\n",
    "YoverLstarfinal = KoverYstarfinal**(alpha/(1-alpha)) * E\n",
    "YoverLstarinitial = KoverYstarinitial**(alpha/(1-alpha)) * E\n",
    "\n",
    "print(KoverYstarfinal, \"= KoverYstarfinal\")\n",
    "print(YoverLstarfinal, \"= YoverLstarfinal\")\n",
    "print(np.log(YoverLstarfinal/YoverLstarinitial), \"= long-run growth effect\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<https://github.com/braddelong/LSS18E101b/blob/master/2018-02-06_Solow_Review.ipynb>\n",
    "\n",
    "<http://datahub.berkeley.edu/user-redirect/interact?account=braddelong&repo=LSS18E101b&branch=master&path=2018-02-06_Solow_Review.ipynb>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}